/*
 * Perfusion.cpp
 *
 *  Created on: 08.03.2011
 *      Author: jagiella
 */

/*
 * main.cpp
 *
 *  Created on: Apr 12, 2010
 *      Author: jagiella
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "Mathematix.h"
#include "VesselGraph.hpp"
#include "EpsIO.h"
//#include "Octree.hpp"

#define MINMAX(  a,  b,  c) \
	( b<a?a: (b>c?c : b))

double rgbformulae(char formulae, double x) {
	double ret = 0;
	//double PI = 3.14159265;
	switch (abs(formulae)) {
	case 0:
		ret = 0.;
		break;
	case 1:
		ret = 0.5;
		break;
	case 2:
		ret = 1.;
		break;
	case 3:
		ret = x;
		break;
	case 4:
		ret = pow(x, 2.);
		break;
	case 5:
		ret = pow(x, 3.);
		break;
	case 6:
		ret = pow(x, 4.);
		break;
	case 7:
		ret = sqrt(x);
		break;
	case 8:
		ret = sqrt(sqrt(x));
		break;
	case 9:
		ret = sin(PI / 2. * x);
		break;
	case 10:
		ret = cos(PI / 2. * x);
		break;
	case 11:
		ret = fabs(x - 0.5);
		break;
	case 12:
		ret = pow(2. * x - 1., 2.);
		break;
	case 13:
		ret = sin(PI * x);
		break;
	case 14:
		ret = fabs(cos(PI * x));
		break;
	case 15:
		ret = sin(2. * PI * x);
		break;
	case 16:
		ret = cos(2. * PI * x);
		break;
	case 17:
		ret = fabs(sin(2. * PI * x));
		break;
	case 18:
		ret = fabs(cos(2. * PI * x));
		break;
	case 19:
		ret = fabs(sin(4. * PI * x));
		break;
	case 20:
		ret = fabs(cos(4. * PI * x));
		break;
	case 21:
		ret = 3. * x;
		break;
	case 22:
		ret = 3. * x - 1.;
		break;
	case 23:
		ret = 3. * x - 2.;
		break;
	case 24:
		ret = fabs(3. * x - 1.);
		break;
	case 25:
		ret = fabs(3. * x - 2.);
		break;
	case 26:
		ret = (3. * x - 1.) / 2.;
		break;
	case 27:
		ret = (3. * x - 2.) / 2.;
		break;
	case 28:
		ret = fabs((3. * x - 1.) / 2.);
		break;
	case 29:
		ret = fabs((3. * x - 2.) / 2.);
		break;
	case 30:
		ret = x / 0.32 - 0.78125;
		break;
	case 31:
		ret = 2. * x - 0.84;
		break;
		//case 32: ret = 4x;1;-2x+1.84;x/0.08-11.5; break;
	case 33:
		ret = fabs(2. * x - 0.5);
		break;
	case 34:
		ret = 2. * x;
		break;
	case 35:
		ret = 2. * x - 0.5;
		break;
	case 36:
		ret = 2. * x - 1.;
		break;
	}

	if (formulae < 0)
		return MINMAX(0., -ret, 1.);
	else
		return MINMAX(0., ret, 1.);
}

/*color_t rgbformulaeMapping(char rformulae, char gformulae, char bformulae,
		double x, double min, double max) {
	color_t color = 0;
	char *formulae[3];
	formulae[0] = &rformulae;
	formulae[1] = &gformulae;
	formulae[2] = &bformulae;

	for (int c = 0; c < 3; c++) {
		//fprintf(stderr, "%lf\n", rgbformulae( *formulae[c], (x-min) / (max-min)));
		((channel_t*) (&color))[c] = (channel_t) 255. * rgbformulae(
				*formulae[c], (x - min) / (max - min));
	}

	return color;
}*/


void printMarkerIntensityMap( char *filename, VesselGraph *vg, float ***markerIntSpace, double voxelSizeX, double voxelSizeY, double voxelSizeZ, bool printAIF, float markerAIF)
{
	// COARSENING
	//double avg = 5; // *60ym = 300ym = 0.3mm
	int X = MAX( 1, DOMAIN_SIZE_X * LATTICE_CONSTANT / voxelSizeX);
	int Y = MAX( 1, DOMAIN_SIZE_Y * LATTICE_CONSTANT / voxelSizeY);
	int Z = MAX( 1, DOMAIN_SIZE_Z * LATTICE_CONSTANT / voxelSizeZ);

	int dx = DOMAIN_SIZE_X/X;
	int dy = DOMAIN_SIZE_Y/Y;
	int dz = DOMAIN_SIZE_Z/Z;

	double ***flowAvg = (double***) malloc(sizeof(double**) * X);
	for (int x = 0; x < X; x++){
		flowAvg[x] = (double**) malloc(sizeof(double*) * Y);
		for (int y = 0; y < Y; y++){
			flowAvg[x][y] = (double*) malloc(sizeof(double) * Z);
			for (int z = 0; z < Z; z++){
				flowAvg[x][y][z] = 0.;
			}
		}
	}

	// Add Marker in Interstitial Space
	float min=0, max=0;
	for (int x = 0; x < DOMAIN_SIZE_X; x++){
		for (int y = 0; y < DOMAIN_SIZE_Y; y++){
			for (int z = 0; z < DOMAIN_SIZE_Z; z++){
				//flowAvg[x/dx][y/dy][z/dz] += markerIntSpace[x][y][z];

				// volume: tissue
				float volume = LATTICE_CONSTANT*LATTICE_CONSTANT*LATTICE_CONSTANT;
				float volumeVessels = 0.;
				float volumeInt = 0.;

				VesselNode *vn = vg->octree->at( x,y,z);
				if( vn!=0){
					// volume of vessels
					for( int v=0; v<vn->countNeighbors; v++)
						volumeVessels += LATTICE_CONSTANT* PI * vn->branches[v]->radius*vn->branches[v]->radius;
					volumeVessels *= 0.5;

					// volume of interstitial space
					volumeInt = volume - volumeVessels;

					flowAvg[x/dx][y/dy][z/dz] +=
						volumeInt    /volume * markerIntSpace[x][y][z] +
						volumeVessels/volume * vn->marker;
				}else{
					flowAvg[x/dx][y/dy][z/dz] += markerIntSpace[x][y][z];
				}
//				if( max<flowAvg[x/dx][y/dy][z/dz])
//					max=flowAvg[x/dx][y/dy][z/dz];

			}
		}
	}

	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
				flowAvg[x][y][z] /= (double)(dx*dy*dz);

	//max = ARTERIAL_MARKER_CONC;

	if( printAIF)
		flowAvg[0][0][0] = markerAIF;

	max = 1.;//(dx*dy*dz);

	// OUTPUT
	std::fstream fs;
	fs.open( filename, std::fstream::out);
	EpsIO::PSwriteHeader(&fs, 0, X, 0, Y);
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
		{
			double xPoly[5] = {(double)x, (double)x+1, (double)x+1, (double)x  ,(double)x};
			double yPoly[5] = {(double)y, (double)y  , (double)y+1, (double)y+1,(double)y};
			char color[512];
			sprintf( color, "%lf %lf %lf setrgbcolor",
					(flowAvg[x][y][z]-min)/(max-min), (flowAvg[x][y][z]-min)/(max-min), (flowAvg[x][y][z]-min)/(max-min));
			EpsIO::PSfillPolygon( &fs,
					xPoly, yPoly, 4,
					0.1, color);
		}
	fs.close();
	for (int x = 0; x < X; x++){
		for (int y = 0; y < Y; y++)
			free(flowAvg[x][y]);
		free(flowAvg[x]);
	}
	free(flowAvg);
}

void printColorMapLegend( char *filename, char *min, char *max)
{
	std::fstream fs;
	fs.open( filename, std::fstream::out);
	EpsIO::PSwriteHeader(&fs, 0, 100, 0,	100);

	for( int i=0; i<100; i++){
		double xPoly[5] = {(double)0, (double)20, (double)20, (double)0  ,(double)0};
		double yPoly[5] = {(double)i, (double)i  , (double)i+1, (double)i+1,(double)i};
		char color[512];
		sprintf( color, "%lf %lf %lf setrgbcolor",
				rgbformulae(  22, i/100.),
				rgbformulae(  13, i/100.),
				rgbformulae( -31, i/100.)
				);

		EpsIO::PSfillPolygon( &fs,
				xPoly, yPoly, 4,
				0.1, color);
	}
	fs<<"0.0 0.0 0.0 setrgbcolor";
	fs<<"/Times-Roman findfont 15 scalefont setfont 20 5 moveto ("<< min <<") show ";
	fs<<"/Times-Roman findfont 15 scalefont setfont 20 85 moveto ("<< max <<") show ";

	fs.close();
}


void printFlowMap( char *filename, VesselGraph *vg, double voxelSizeX, double voxelSizeY, double voxelSizeZ)
{
	// COARSENING
	//double avg = 5; // *60ym = 300ym = 0.3mm

	// domain size (in voxels)
	int X = MAX( 1, DOMAIN_SIZE_X * LATTICE_CONSTANT / voxelSizeX);
	int Y = MAX( 1, DOMAIN_SIZE_Y * LATTICE_CONSTANT / voxelSizeY);
	int Z = MAX( 1, DOMAIN_SIZE_Z * LATTICE_CONSTANT / voxelSizeZ);

	// voxel size (in lattice sites)
	int dx = DOMAIN_SIZE_X/X;
	int dy = DOMAIN_SIZE_Y/Y;
	int dz = DOMAIN_SIZE_Z/Z;

	// INIT AVERAGE DATA STRUCTURE
	double ***flowAvg = (double***) malloc(sizeof(double**) * X);
	for (int x = 0; x < X; x++){
		flowAvg[x] = (double**) malloc(sizeof(double*) * Y);
		for (int y = 0; y < Y; y++){
			flowAvg[x][y] = (double*) malloc(sizeof(double) * Z);
			for (int z = 0; z < Z; z++){
				flowAvg[x][y][z] = 0.;
			}
		}
	}
	double ***flowAvgOut = (double***) malloc(sizeof(double**) * X);
	for (int x = 0; x < X; x++){
		flowAvgOut[x] = (double**) malloc(sizeof(double*) * Y);
		for (int y = 0; y < Y; y++){
			flowAvgOut[x][y] = (double*) malloc(sizeof(double) * Z);
			for (int z = 0; z < Z; z++){
				flowAvgOut[x][y][z] = 0.;
			}
		}
	}

	// SUM UP VALUES PER VOXEL
	float min=0, max=0;
	for (int x = 0; x < DOMAIN_SIZE_X; x++){
		for (int y = 0; y < DOMAIN_SIZE_Y; y++){
			for (int z = 0; z < DOMAIN_SIZE_Z; z++){

				// actual lattice site
				VesselNode *vn = vg->octree->at( x,y,z);

				if(vn!=0){
				// out flow into neighboring voxel?
				for( int v=0; v<vn->countNeighbors; v++)
				{
					//double outFlow = 0.;
					//double inFlow = 0.;
					// out flow?
					if( vn->pressure >= vn->neighbors[v]->pressure){
						// yes!

						// neighbor in neighboring voxel?
						if( vn->neighbors[v]->position[0]/dx != x/dx ||
							vn->neighbors[v]->position[1]/dy != y/dy ||
							vn->neighbors[v]->position[2]/dz != z/dz){
							// yes!
							flowAvgOut[x/dx][y/dy][z/dz] += fabs(vn->branches[v]->flow);
							//outFlow += fabs(vn->branches[v]->flow);
						}
					}
					// in flow!
					else{
						// neighbor in neighboring voxel?
						if( vn->neighbors[v]->position[0]/dx != x/dx ||
							vn->neighbors[v]->position[1]/dy != y/dy ||
							vn->neighbors[v]->position[2]/dz != z/dz){
							// yes!
							flowAvg[x/dx][y/dy][z/dz] += fabs(vn->branches[v]->flow);
						}
					}
					//flowAvg[x/dx][y/dy][z/dz] += MAX(inFlow,outFlow);
				}
				}
			}
		}
	}

	// NORMALIZE
	/*for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
				flowAvg[x][y][z] /= (double)(dx*dy*dz);
	 */

	// MIN & MAX
	min=MAX(flowAvg[0][0][0], flowAvgOut[0][0][0]);
	max=MAX(flowAvg[0][0][0], flowAvgOut[0][0][0]);
	for (int x = 0; x < X; x++)
			for (int y = 0; y < Y; y++)
				for (int z = 0; z < 1; z++){
					flowAvg[x][y][z] = MAX(flowAvg[x][y][z], flowAvgOut[x][y][z]);
					if( min>flowAvg[x][y][z])
						min=flowAvg[x][y][z];

					if( max<flowAvg[x][y][z])
						max=flowAvg[x][y][z];
				}

	// OUTPUT
	std::fstream fs;
	fs.open( filename, std::fstream::out);
	EpsIO::PSwriteHeader(&fs, 0, X, 0, Y);
	if(min == max)
		min = 0;
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
		{
			double xPoly[5] = {(double)x, (double)x+1, (double)x+1, (double)x  ,(double)x};
			double yPoly[5] = {(double)y, (double)y  , (double)y+1, (double)y+1,(double)y};
			char color[512];
			sprintf( color, "%lf %lf %lf setrgbcolor",
					//(flowAvg[x][y][z]-min)/(max-min),
					//(flowAvg[x][y][z]-min)/(max-min),
					//(flowAvg[x][y][z]-min)/(max-min)
					rgbformulae(  22, (flowAvg[x][y][z]-min)/(max-min)),
					rgbformulae(  13, (flowAvg[x][y][z]-min)/(max-min)),
					rgbformulae( -31, (flowAvg[x][y][z]-min)/(max-min))
					);
			EpsIO::PSfillPolygon( &fs,
					xPoly, yPoly, 4,
					0.1, color);
		}
	fs.close();

	// LEGEND
	char legend_filename[512];
	char min_string[512];
	char max_string[512];
	sprintf(legend_filename, "%s.legend.eps", filename);
	sprintf(min_string, "%lf", min);
	sprintf(max_string, "%lf", max);
	printColorMapLegend( legend_filename, min_string, max_string);

	// DATA
	char data_filename[512];
	sprintf(data_filename, "%s.data.dat", filename);
	FILE * fp = fopen(data_filename,"w+");
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
				fprintf(fp, "%i %i %e\n", x, y, flowAvg[x][y][z]);
	fclose(fp);

	for (int x = 0; x < X; x++){
		for (int y = 0; y < Y; y++)
			free(flowAvg[x][y]);
		free(flowAvg[x]);
	}
	free(flowAvg);
	for (int x = 0; x < X; x++){
		for (int y = 0; y < Y; y++)
			free(flowAvgOut[x][y]);
		free(flowAvgOut[x]);
	}
	free(flowAvgOut);

}


void printPermeabilityMap( char *filename, VesselGraph *vg, double voxelSizeX, double voxelSizeY, double voxelSizeZ)
{
	// COARSENING
//	double avg = 5; // *60ym = 300ym = 0.3mm

	// domain size (in voxels)
	int X = MAX( 1, DOMAIN_SIZE_X * LATTICE_CONSTANT / voxelSizeX);
	int Y = MAX( 1, DOMAIN_SIZE_Y * LATTICE_CONSTANT / voxelSizeY);
	int Z = MAX( 1, DOMAIN_SIZE_Z * LATTICE_CONSTANT / voxelSizeZ);

	// voxel size (in lattice sites)
	int dx = DOMAIN_SIZE_X/X;
	int dy = DOMAIN_SIZE_Y/Y;
	int dz = DOMAIN_SIZE_Z/Z;

	// INIT AVERAGE DATA STRUCTURE
	double ***flowAvg = (double***) malloc(sizeof(double**) * X);
	for (int x = 0; x < X; x++){
		flowAvg[x] = (double**) malloc(sizeof(double*) * Y);
		for (int y = 0; y < Y; y++){
			flowAvg[x][y] = (double*) malloc(sizeof(double) * Z);
			for (int z = 0; z < Z; z++){
				flowAvg[x][y][z] = 0.;
			}
		}
	}

	// SUM UP VALUES PER VOXEL
	float min=0, max=0;
	for (int x = 0; x < DOMAIN_SIZE_X; x++){
		for (int y = 0; y < DOMAIN_SIZE_Y; y++){
			for (int z = 0; z < DOMAIN_SIZE_Z; z++){

				// actual lattice site
				VesselNode *vn = vg->octree->at( x,y,z);

				if(vn!=0){
					flowAvg[x/dx][y/dy][z/dz] += vg->getVascularPermeabilitySurfaceProduct(vn);
				}
			}
		}
	}

	// NORMALIZE
	/*for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
				flowAvg[x][y][z] /= (double)(dx*dy*dz);
	 */

	// MIN & MAX
	min=flowAvg[0][0][0];
	max=flowAvg[0][0][0];
	for (int x = 0; x < X; x++)
			for (int y = 0; y < Y; y++)
				for (int z = 0; z < 1; z++){
					if( min>flowAvg[x][y][z])
						min=flowAvg[x][y][z];

					if( max<flowAvg[x][y][z])
						max=flowAvg[x][y][z];
				}

	// OUTPUT
	std::fstream fs;
	fs.open( filename, std::fstream::out);
	EpsIO::PSwriteHeader(&fs, 0, X, 0, Y);
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
		{
			double xPoly[5] = {(double)x, (double)x+1, (double)x+1, (double)x  ,(double)x};
			double yPoly[5] = {(double)y, (double)y  , (double)y+1, (double)y+1,(double)y};
			char color[512];
			sprintf( color, "%lf %lf %lf setrgbcolor",
					//(flowAvg[x][y][z]-min)/(max-min),
					//(flowAvg[x][y][z]-min)/(max-min),
					//(flowAvg[x][y][z]-min)/(max-min)
					rgbformulae(  22, (flowAvg[x][y][z]-min)/(max-min)),
					rgbformulae(  13, (flowAvg[x][y][z]-min)/(max-min)),
					rgbformulae( -31, (flowAvg[x][y][z]-min)/(max-min))
					);
			EpsIO::PSfillPolygon( &fs,
					xPoly, yPoly, 4,
					0.1, color);
		}
	fs.close();

	// LEGEND
	char legend_filename[512];
	char min_string[512];
	char max_string[512];
	sprintf(legend_filename, "%s.legend.eps", filename);
	sprintf(min_string, "%lf", min);
	sprintf(max_string, "%lf", max);
	printColorMapLegend( legend_filename, min_string, max_string);

	// DATA
	char data_filename[512];
	sprintf(data_filename, "%s.data.dat", filename);
	FILE * fp = fopen(data_filename,"w+");
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
				fprintf(fp, "%i %i %e\n", x, y, flowAvg[x][y][z]);
	fclose(fp);

	for (int x = 0; x < X; x++){
		for (int y = 0; y < Y; y++)
			free(flowAvg[x][y]);
		free(flowAvg[x]);
	}
	free(flowAvg);
}


void printBloodVolumeFractionMap( char *filename, VesselGraph *vg, double voxelSizeX, double voxelSizeY, double voxelSizeZ)
{
	// COARSENING
//	double avg = 5; // *60ym = 300ym = 0.3mm

	// domain size (in voxels)
	int X = MAX( 1, DOMAIN_SIZE_X * LATTICE_CONSTANT / voxelSizeX);
	int Y = MAX( 1, DOMAIN_SIZE_Y * LATTICE_CONSTANT / voxelSizeY);
	int Z = MAX( 1, DOMAIN_SIZE_Z * LATTICE_CONSTANT / voxelSizeZ);

	// voxel size (in lattice sites)
	int dx = DOMAIN_SIZE_X/X;
	int dy = DOMAIN_SIZE_Y/Y;
	int dz = DOMAIN_SIZE_Z/Z;

	// INIT AVERAGE DATA STRUCTURE
	double ***flowAvg = (double***) malloc(sizeof(double**) * X);
	for (int x = 0; x < X; x++){
		flowAvg[x] = (double**) malloc(sizeof(double*) * Y);
		for (int y = 0; y < Y; y++){
			flowAvg[x][y] = (double*) malloc(sizeof(double) * Z);
			for (int z = 0; z < Z; z++){
				flowAvg[x][y][z] = 0.;
			}
		}
	}

	// SUM UP VALUES PER VOXEL
	float min=0, max=0;
	for (int x = 0; x < DOMAIN_SIZE_X; x++){
		for (int y = 0; y < DOMAIN_SIZE_Y; y++){
			for (int z = 0; z < DOMAIN_SIZE_Z; z++){

				// actual lattice site
				flowAvg[x/dx][y/dy][z/dz] += (double)vg->getVascularVolume( x,y,z)/(LATTICE_CONSTANT*LATTICE_CONSTANT*LATTICE_CONSTANT);
			}
		}
	}

	// NORMALIZE
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
				flowAvg[x][y][z] /= (double)(dx*dy*dz);


	// MIN & MAX
	min=flowAvg[0][0][0];
	max=flowAvg[0][0][0];
	for (int x = 0; x < X; x++)
			for (int y = 0; y < Y; y++)
				for (int z = 0; z < 1; z++){
					if( min>flowAvg[x][y][z])
						min=flowAvg[x][y][z];

					if( max<flowAvg[x][y][z])
						max=flowAvg[x][y][z];
				}

	// OUTPUT
	std::fstream fs;
	fs.open( filename, std::fstream::out);
	EpsIO::PSwriteHeader(&fs, 0, X, 0, Y);
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
		{
			double xPoly[5] = {(double)x, (double)x+1, (double)x+1, (double)x  ,(double)x};
			double yPoly[5] = {(double)y, (double)y  , (double)y+1, (double)y+1,(double)y};
			char color[512];
			sprintf( color, "%lf %lf %lf setrgbcolor",
					//(flowAvg[x][y][z]-min)/(max-min),
					//(flowAvg[x][y][z]-min)/(max-min),
					//(flowAvg[x][y][z]-min)/(max-min)
					rgbformulae(  22, (flowAvg[x][y][z]-min)/(max-min)),
					rgbformulae(  13, (flowAvg[x][y][z]-min)/(max-min)),
					rgbformulae( -31, (flowAvg[x][y][z]-min)/(max-min))
					);
			EpsIO::PSfillPolygon( &fs,
					xPoly, yPoly, 4,
					0.1, color);
		}
	fs.close();

	char legend_filename[512];
	char min_string[512];
	char max_string[512];
	sprintf(legend_filename, "%s.legend.eps", filename);
	sprintf(min_string, "%f", min);
	sprintf(max_string, "%f", max);
	printColorMapLegend( legend_filename, min_string, max_string);

	// DATA
	char data_filename[512];
	sprintf(data_filename, "%s.data.dat", filename);
	FILE * fp = fopen(data_filename,"w+");
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
				fprintf(fp, "%i %i %e\n", x, y, flowAvg[x][y][z]);
	fclose(fp);

	for (int x = 0; x < X; x++){
		for (int y = 0; y < Y; y++)
			free(flowAvg[x][y]);
		free(flowAvg[x]);
	}
	free(flowAvg);
}


void printInterstitialSpaceVolumeFractionMap( char *filename, VesselGraph *vg, double voxelSizeX, double voxelSizeY, double voxelSizeZ)
{
	// COARSENING
//	double avg = 5; // *60ym = 300ym = 0.3mm

	// domain size (in voxels)
	int X = MAX( 1, DOMAIN_SIZE_X * LATTICE_CONSTANT / voxelSizeX);
	int Y = MAX( 1, DOMAIN_SIZE_Y * LATTICE_CONSTANT / voxelSizeY);
	int Z = MAX( 1, DOMAIN_SIZE_Z * LATTICE_CONSTANT / voxelSizeZ);

	// voxel size (in lattice sites)
	int dx = DOMAIN_SIZE_X/X;
	int dy = DOMAIN_SIZE_Y/Y;
	int dz = DOMAIN_SIZE_Z/Z;

	// INIT AVERAGE DATA STRUCTURE
	double ***flowAvg = (double***) malloc(sizeof(double**) * X);
	for (int x = 0; x < X; x++){
		flowAvg[x] = (double**) malloc(sizeof(double*) * Y);
		for (int y = 0; y < Y; y++){
			flowAvg[x][y] = (double*) malloc(sizeof(double) * Z);
			for (int z = 0; z < Z; z++){
				flowAvg[x][y][z] = 0.;
			}
		}
	}

	// SUM UP VALUES PER VOXEL
	float min=0, max=0;
	for (int x = 0; x < DOMAIN_SIZE_X; x++){
		for (int y = 0; y < DOMAIN_SIZE_Y; y++){
			for (int z = 0; z < DOMAIN_SIZE_Z; z++){

				// actual lattice site
				flowAvg[x/dx][y/dy][z/dz] += (double)vg->getExtraVascularVolume( x,y,z)/(LATTICE_CONSTANT*LATTICE_CONSTANT*LATTICE_CONSTANT);
			}
		}
	}

	// NORMALIZE
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
				flowAvg[x][y][z] /= (double)(dx*dy*dz);


	// MIN & MAX
	min=flowAvg[0][0][0];
	max=flowAvg[0][0][0];
	for (int x = 0; x < X; x++)
			for (int y = 0; y < Y; y++)
				for (int z = 0; z < 1; z++){
					if( min>flowAvg[x][y][z])
						min=flowAvg[x][y][z];

					if( max<flowAvg[x][y][z])
						max=flowAvg[x][y][z];
				}

	// OUTPUT
	std::fstream fs;
	fs.open( filename, std::fstream::out);
	EpsIO::PSwriteHeader(&fs, 0, X, 0, Y);
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
		{
			double xPoly[5] = {(double)x, (double)x+1, (double)x+1, (double)x  ,(double)x};
			double yPoly[5] = {(double)y, (double)y  , (double)y+1, (double)y+1,(double)y};
			char color[512];
			sprintf( color, "%lf %lf %lf setrgbcolor",
					//(flowAvg[x][y][z]-min)/(max-min),
					//(flowAvg[x][y][z]-min)/(max-min),
					//(flowAvg[x][y][z]-min)/(max-min)
					rgbformulae(  22, (flowAvg[x][y][z]-min)/(max-min)),
					rgbformulae(  13, (flowAvg[x][y][z]-min)/(max-min)),
					rgbformulae( -31, (flowAvg[x][y][z]-min)/(max-min))
					);
			EpsIO::PSfillPolygon( &fs,
					xPoly, yPoly, 4,
					0.1, color);
		}
	fs.close();

	char legend_filename[512];
	char min_string[512];
	char max_string[512];
	sprintf(legend_filename, "%s.legend.eps", filename);
	sprintf(min_string, "%f", min);
	sprintf(max_string, "%f", max);
	printColorMapLegend( legend_filename, min_string, max_string);

	// DATA
	char data_filename[512];
	sprintf(data_filename, "%s.data.dat", filename);
	FILE * fp = fopen(data_filename,"w+");
	for (int x = 0; x < X; x++)
		for (int y = 0; y < Y; y++)
			for (int z = 0; z < 1; z++)
				fprintf(fp, "%i %i %e\n", x, y, flowAvg[x][y][z]);
	fclose(fp);

	for (int x = 0; x < X; x++){
		for (int y = 0; y < Y; y++)
			free(flowAvg[x][y]);
		free(flowAvg[x]);
	}
	free(flowAvg);
}


VesselNode *getClosestVesselNode(VesselGraph *vg, float *pos, float &sqrDist) {
	if (vg->countVesselNodes == 0)
		return 0;

	float dist;
	sqrDist = pow(pos[0] - vg->vesselNodes[0]->position[0], 2) + pow(pos[1]
			- vg->vesselNodes[0]->position[1], 2) + pow(pos[2]
			- vg->vesselNodes[0]->position[2], 2);
	VesselNode *closestVN = vg->vesselNodes[0];
	for (int v = 1; v < vg->countVesselNodes; v++) {
		dist = pow(pos[0] - vg->vesselNodes[v]->position[0], 2) + pow(pos[1]
				- vg->vesselNodes[v]->position[1], 2) + pow(pos[2]
				- vg->vesselNodes[v]->position[2], 2);
		if (sqrDist > dist) {
			sqrDist = dist;
			closestVN = vg->vesselNodes[v];
		}
	}

	return closestVN;
}

void setInitialVesselGraph(VesselGraph *vesselGraph) {


	// PLACE ROOT NODES
	float position[3];
	float dist = 0.;
	float pressure = MAX_PRESSURE;//MIN_PRESSURE;

	//fprintf( stderr, "Add Root Node (%i, %i)\n", (int)DOMAIN_SIZE/2,(int)DOMAIN_SIZE/2);
#if DIMENSIONS == 1
	 //VesselNode *last = new VesselNode(DOMAIN_SIZE_X/2,0,0);
	 VesselNode *last = new VesselNode(0,0,0);
#elif DIMENSIONS == 2
	 VesselNode *last = new VesselNode(DOMAIN_SIZE_X/2,DOMAIN_SIZE_Y/2,0);
#else
	 VesselNode *last = new VesselNode(DOMAIN_SIZE_X/2,DOMAIN_SIZE_Y/2,0);
#endif
	 vesselGraph->addVesselNode( last);
	 last->pressure = pressure; // kPa
	 last->setType( ROOT);
	fprintf(stderr, "Added Root Node (%.0f, %.0f, %.0f) -> pressure = %lf kPa\n",
			last->position[0],last->position[1],last->position[2],
			last->pressure);


	for (int i = 0; i < 1000; i++) {
		int dx[3] = { 0, 0, 0 };

		// RANDOM POSITION AT DOMAIN BORDER
		switch (i % (2 * DIMENSIONS)) {
		case 0:
			// YZ-AXIS
			position[0] = 0;
			position[1] = BRANCH_LENGTH * floor((DOMAIN_SIZE_Y - 1)	/ BRANCH_LENGTH * RAND01);
			position[2] = BRANCH_LENGTH * floor((DOMAIN_SIZE_Z - 1)	/ BRANCH_LENGTH * RAND01);
			dx[0] = 1;
			break;
		case 1:
			// YZ-AXIS
			position[0] = DOMAIN_SIZE_X - 1;
			position[1] = BRANCH_LENGTH * floor((DOMAIN_SIZE_Y - 1)	/ BRANCH_LENGTH * RAND01);
			position[2] = BRANCH_LENGTH * floor((DOMAIN_SIZE_Z - 1)	/ BRANCH_LENGTH * RAND01);
			dx[0] = -1;
			break;
		case 2:
			// XZ-AXIS
			position[0] = BRANCH_LENGTH * floor((DOMAIN_SIZE_X - 1)	/ BRANCH_LENGTH * RAND01);
			position[1] = 0;
			position[2] = BRANCH_LENGTH * floor((DOMAIN_SIZE_Z - 1)	/ BRANCH_LENGTH * RAND01);
			dx[1] = 1;
			break;
		case 3:
			// XZ-AXIS
			position[0] = BRANCH_LENGTH * floor((DOMAIN_SIZE_X - 1)
					/ BRANCH_LENGTH * RAND01);
			position[1] = DOMAIN_SIZE_Y - 1;
			position[2] = BRANCH_LENGTH * floor((DOMAIN_SIZE_Z - 1)
					/ BRANCH_LENGTH * RAND01);
			dx[1] = -1;
			break;
			//#if DIMENSIONS==3
		case 4:
			// XY-AXIS
			position[0] = BRANCH_LENGTH * floor((DOMAIN_SIZE_X - 1)
					/ BRANCH_LENGTH * RAND01);
			position[1] = BRANCH_LENGTH * floor((DOMAIN_SIZE_Y - 1)
					/ BRANCH_LENGTH * RAND01);
			position[2] = 0;
			dx[2] = 1;
			break;
		case 5:
			// XY-AXIS
			position[0] = BRANCH_LENGTH * floor((DOMAIN_SIZE_X - 1)
					/ BRANCH_LENGTH * RAND01);
			position[1] = BRANCH_LENGTH * floor((DOMAIN_SIZE_Y - 1)
					/ BRANCH_LENGTH * RAND01);
			position[2] = DOMAIN_SIZE_Z - 1;
			dx[2] = -1;
			break;
			//#endif

		}
		VesselNode *closest = getClosestVesselNode(vesselGraph, position, dist);
		//if (closest != 0)
//			pressure = (closest->pressure == MIN_PRESSURE ? MAX_PRESSURE : MIN_PRESSURE);

		if (vesselGraph->countVesselNodes == 0 || dist >= ROOT_DISTANCE*ROOT_DISTANCE) {
			if (closest != 0)
				fprintf(stderr, "Dist: %lf, Type: %i, Pressure: %lf\n", dist, closest->getType(), closest->pressure);
			VesselNode *last = new VesselNode(position[0], position[1],
					position[2]);
			vesselGraph->addVesselNode(last);
			last->pressure = MIN_PRESSURE;//MAX_PRESSURE;//pressure; // kPa
			//last->pressure = pressure; // kPa
			last->setType(ROOT);
			fprintf(stderr, "Added Root Node (%i, %i, %i) -> pressure = %lf kPa\n",
					(int) position[0], (int) position[1], (int) position[2],
					last->pressure);

			VesselNode *next;
			for (int i = 0; i < BRANCH_LENGTH; i++) {
				next = new VesselNode(last->position[0] + dx[0],
						last->position[1] + dx[1], last->position[2] + dx[2]);
				next->setType(VESSEL);
				vesselGraph->addVesselNode(next);
				vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(last, next));
				last = next;
			}
			next->setType(TIP);

			pressure = (pressure == MIN_PRESSURE ? MAX_PRESSURE : MIN_PRESSURE);
		}
	}

	// GROW VESSEL GRAPH STOCHASTICALLY
	bool changed;
	do {
		changed = false;
		int countPreviousVesselNodes = vesselGraph->countVesselNodes;
		for (int v = 0; v < countPreviousVesselNodes; v++)
			//if (RAND01 < 0.5)
			{

				// BRANCHING
				if ((vesselGraph->vesselNodes[v]->getType() == TIP /*|| vesselGraph->vesselNodes[v]->getType() == ROOT*/)
						&& vesselGraph->vesselNodes[v]->countNeighbors <= 1) {
					//fprintf( stderr, "BRANCH?\n");
					int x = vesselGraph->vesselNodes[v]->position[0];
					int y = vesselGraph->vesselNodes[v]->position[1];
					int z = vesselGraph->vesselNodes[v]->position[2];
					if (x >= BRANCH_LENGTH && vesselGraph->octree->at(x
							- BRANCH_LENGTH, y, z) == 0 && x < DOMAIN_SIZE_X
							- BRANCH_LENGTH && vesselGraph->octree->at(x
							+ BRANCH_LENGTH, y, z) == 0) {
						//fprintf( stderr, "BRANCH\n");
						changed = true;
						if (vesselGraph->vesselNodes[v]->getType() != ROOT)
							vesselGraph->vesselNodes[v]->setType(VESSEL);
						VesselNode *last, *next;

						int dx[3] = { 1, 0, 0 };
						last = vesselGraph->vesselNodes[v];
						for (int i = 0; i < BRANCH_LENGTH; i++) {
							next = new VesselNode(
									vesselGraph->vesselNodes[v]->position[0]
											+ (i + 1) * dx[0],
									vesselGraph->vesselNodes[v]->position[1]
											+ (i + 1) * dx[1],
									vesselGraph->vesselNodes[v]->position[2]
											+ (i + 1) * dx[2]);
							next->setType(VESSEL);
							vesselGraph->addVesselNode(next);
							vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment( last, next));
							last = next;
						}
						next->setType(TIP);

						dx[0] = -1;
						last = vesselGraph->vesselNodes[v];
						for (int i = 0; i < BRANCH_LENGTH; i++) {
							next = new VesselNode(
									vesselGraph->vesselNodes[v]->position[0]
											+ (i + 1) * dx[0],
									vesselGraph->vesselNodes[v]->position[1]
											+ (i + 1) * dx[1],
									vesselGraph->vesselNodes[v]->position[2]
											+ (i + 1) * dx[2]);
							next->setType(VESSEL);
							vesselGraph->addVesselNode(next);
							vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
									last, next));
							last = next;
						}
						next->setType(TIP);
					} else if (y >= BRANCH_LENGTH
							&& vesselGraph->octree->at(x, y - BRANCH_LENGTH,z) == 0
							&& y < DOMAIN_SIZE_Y - BRANCH_LENGTH
							&& vesselGraph->octree->at(x, y + BRANCH_LENGTH, z) == 0) {
						//fprintf( stderr, "BRANCH\n");
						changed = true;
						if (vesselGraph->vesselNodes[v]->getType() != ROOT)
							vesselGraph->vesselNodes[v]->setType(VESSEL);
						VesselNode *last, *next;

						int dx[3] = { 0, 1, 0 };
						last = vesselGraph->vesselNodes[v];
						for (int i = 0; i < BRANCH_LENGTH; i++) {
							next = new VesselNode(
									vesselGraph->vesselNodes[v]->position[0]
											+ (i + 1) * dx[0],
									vesselGraph->vesselNodes[v]->position[1]
											+ (i + 1) * dx[1],
									vesselGraph->vesselNodes[v]->position[2]
											+ (i + 1) * dx[2]);
							next->setType(VESSEL);
							vesselGraph->addVesselNode(next);
							vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
									last, next));
							last = next;
						}
						next->setType(TIP);

						dx[1] = -1;
						last = vesselGraph->vesselNodes[v];
						for (int i = 0; i < BRANCH_LENGTH; i++) {
							next = new VesselNode(
									vesselGraph->vesselNodes[v]->position[0]
											+ (i + 1) * dx[0],
									vesselGraph->vesselNodes[v]->position[1]
											+ (i + 1) * dx[1],
									vesselGraph->vesselNodes[v]->position[2]
											+ (i + 1) * dx[2]);
							next->setType(VESSEL);
							vesselGraph->addVesselNode(next);
							vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
									last, next));
							last = next;
						}
						next->setType(TIP);
					} else if (z >= BRANCH_LENGTH
							&& vesselGraph->octree->at(x, y, z
									- BRANCH_LENGTH) == 0 && z
							< DOMAIN_SIZE_Z - BRANCH_LENGTH
							&& vesselGraph->octree->at(x, y, z
									+ BRANCH_LENGTH) == 0) {
						//fprintf( stderr, "BRANCH\n");
						changed = true;
						if (vesselGraph->vesselNodes[v]->getType() != ROOT)
							vesselGraph->vesselNodes[v]->setType(VESSEL);
						VesselNode *last, *next;

						int dx[3] = { 0, 0, 1 };
						last = vesselGraph->vesselNodes[v];
						for (int i = 0; i < BRANCH_LENGTH; i++) {
							next = new VesselNode(
									vesselGraph->vesselNodes[v]->position[0]
											+ (i + 1) * dx[0],
									vesselGraph->vesselNodes[v]->position[1]
											+ (i + 1) * dx[1],
									vesselGraph->vesselNodes[v]->position[2]
											+ (i + 1) * dx[2]);
							next->setType(VESSEL);
							vesselGraph->addVesselNode(next);
							vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
									last, next));
							last = next;
						}
						next->setType(TIP);

						dx[2] = -1;
						last = vesselGraph->vesselNodes[v];
						for (int i = 0; i < BRANCH_LENGTH; i++) {
							next = new VesselNode(
									vesselGraph->vesselNodes[v]->position[0]
											+ (i + 1) * dx[0],
									vesselGraph->vesselNodes[v]->position[1]
											+ (i + 1) * dx[1],
									vesselGraph->vesselNodes[v]->position[2]
											+ (i + 1) * dx[2]);
							next->setType(VESSEL);
							vesselGraph->addVesselNode(next);
							vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
									last, next));
							last = next;
						}
						next->setType(TIP);
					}
				}

				// SINGLE SPROUT
				if ((vesselGraph->vesselNodes[v]->getType() == ROOT	&& vesselGraph->vesselNodes[v]->countNeighbors == 0)
						|| (vesselGraph->vesselNodes[v]->getType() == TIP && vesselGraph->vesselNodes[v]->countNeighbors < 2)) {
					int countCandidates = 0;
					int candidateX[8];
					int candidateY[8];
					int candidateZ[8];
					int x = vesselGraph->vesselNodes[v]->position[0];
					int y = vesselGraph->vesselNodes[v]->position[1];
					int z = vesselGraph->vesselNodes[v]->position[2];
					//float maxVEGF = vegf[x][y];

					for (int xx = x - BRANCH_LENGTH; xx <= x
							+ BRANCH_LENGTH; xx += 2 * BRANCH_LENGTH)
						if (xx >= 0 && xx < DOMAIN_SIZE_X
						//&& (y<BRANCH_LENGTH || vesselGraph->octree->at(xx,y-BRANCH_LENGTH,0) == 0)
								&& (vesselGraph->octree->at(xx, y, z) == 0)
						//&& (y>=DOMAIN_SIZE-BRANCH_LENGTH || vesselGraph->octree->at(xx,y+BRANCH_LENGTH,0) == 0)
						) {
							//maxVEGF = vegf[xx][y];
							candidateX[countCandidates] = xx;
							candidateY[countCandidates] = y;
							candidateZ[countCandidates] = z;
							countCandidates++;
						}

					for (int yy = y - BRANCH_LENGTH; yy <= y
							+ BRANCH_LENGTH; yy += 2 * BRANCH_LENGTH)
						if (yy >= 0 && yy < DOMAIN_SIZE_Y
						//&& (x==0 || vesselGraph->octree->at(x-BRANCH_LENGTH,yy,0) == 0)
								&& vesselGraph->octree->at(x, yy, z) == 0
						//&& (x>=DOMAIN_SIZE-BRANCH_LENGTH || vesselGraph->octree->at(x+BRANCH_LENGTH,yy,0) == 0)
						) {
							//maxVEGF = vegf[x][yy];
							candidateX[countCandidates] = x;
							candidateY[countCandidates] = yy;
							candidateZ[countCandidates] = z;
							countCandidates++;
						}

					for (int zz = z - BRANCH_LENGTH; zz <= z
							+ BRANCH_LENGTH; zz += 2 * BRANCH_LENGTH)
						if (zz >= 0 && zz < DOMAIN_SIZE_Z
						//&& (x==0 || vesselGraph->octree->at(x-BRANCH_LENGTH,yy,0) == 0)
								&& vesselGraph->octree->at(x, y, zz) == 0
						//&& (x>=DOMAIN_SIZE-BRANCH_LENGTH || vesselGraph->octree->at(x+BRANCH_LENGTH,yy,0) == 0)
						) {
							//maxVEGF = vegf[x][yy];
							candidateX[countCandidates] = x;
							candidateY[countCandidates] = y;
							candidateZ[countCandidates] = zz;
							countCandidates++;
						}
					//fprintf( stderr, "Sprout! maxVEGF=%e\n", maxVEGF);

					VesselNode *last, *next;
					if (countCandidates > 0) {
						last = vesselGraph->vesselNodes[v];
						if (vesselGraph->vesselNodes[v]->getType() != ROOT)
							vesselGraph->vesselNodes[v]->setType(VESSEL);
						int which = RAND01 * countCandidates;
						int dx[3] = { (candidateX[which]
								- vesselGraph->vesselNodes[v]->position[0])
								/ BRANCH_LENGTH, (candidateY[which]
								- vesselGraph->vesselNodes[v]->position[1])
								/ BRANCH_LENGTH, (candidateZ[which]
								- vesselGraph->vesselNodes[v]->position[2])
								/ BRANCH_LENGTH };
						for (int i = 0; i < BRANCH_LENGTH; i++) {
							next = new VesselNode(
									vesselGraph->vesselNodes[v]->position[0]
											+ (i + 1) * dx[0],
									vesselGraph->vesselNodes[v]->position[1]
											+ (i + 1) * dx[1],
									vesselGraph->vesselNodes[v]->position[2]
											+ (i + 1) * dx[2]);
							next->setType(VESSEL);
							vesselGraph->addVesselNode(next);
							vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
									last, next));
							last = next;
						}
						next->setType(TIP);
						changed = true;
					}

				}
			}
	} while (changed);
}

void setInitialVesselGraphOLD(VesselGraph *vesselGraph) {
	VesselNode *last;
	VesselNode *next = new VesselNode(0, 0, 0);
	vesselGraph->addVesselNode(next);
	next->pressure = 12.; // kPa
	next->setType(ROOT);
	int position[3] = { 0, 0, 0 };
	int max[3] = { DOMAIN_SIZE_X, DOMAIN_SIZE_Y, DOMAIN_SIZE_Z };
	for (int i = 0; i < DOMAIN_SIZE_X + DOMAIN_SIZE_Y + DOMAIN_SIZE_Z - 2; i++) {

		// chose direction
		//int dx[3] = {0,0,0};
		//dx[(int)(RAND01*3.)] ++;
		int dxi;
		do {
			dxi = (int) (RAND01 * (float) DIMENSIONS);
			fprintf(stderr, "i: %i, dxi: %i, position[dxi]: %i\n", i, dxi,
					position[dxi]);
		} while (position[dxi] == max[dxi] - 1);

		position[dxi]++;
		//fprintf( stderr, "(%i, %i, %i)\n", position[0], position[1], position[2]);

		// new vessel node
		last = next;
		next = new VesselNode(position[0], position[1], position[2]);
		next->setType(VESSEL);
		vesselGraph->addVesselNode(next);
		vesselGraph->addVesselSegment(new VesselSegment(last, next));
		vesselGraph->vesselSegments[vesselGraph->countVesselSegments - 1]->radius
				= 25;

		position[dxi]++;
		last = next;
		next = new VesselNode(position[0], position[1], position[2]);
		next->setType(VESSEL);
		vesselGraph->addVesselNode(next);
		vesselGraph->addVesselSegment(new VesselSegment(last, next));
		vesselGraph->vesselSegments[vesselGraph->countVesselSegments - 1]->radius
				= 25;
		i++;

		//fprintf( stdout, "%i %i %i %i %i %i\n",
		//		last->position[0], last->position[1], last->position[2],
		//		next->position[0]-last->position[0], next->position[1]-last->position[1], next->position[2]-last->position[2]);
	}
	next->pressure = 12;
	next->setType(ROOT);

	next = new VesselNode(0, 100, 0);
	vesselGraph->addVesselNode(next);
	next->pressure = 0.; // kPa
	next->setType(ROOT);
	position[0] = 0;
	position[1] = 100;
	position[2] = 0;
	for (int i = 0; i < DOMAIN_SIZE_X + DOMAIN_SIZE_Y + DOMAIN_SIZE_Z - 2; i++) {

		// chose direction
		//int dx[3] = {0,0,0};
		//dx[(int)(RAND01*3.)] ++;
		int dxi;
		do {
			dxi = (int) (RAND01 * (float) DIMENSIONS);
			fprintf(stderr, "i: %i, dxi: %i, position[dxi]: %i\n", i, dxi,
					position[dxi]);
		} while ((dxi == 0 && position[dxi] == max[dxi] - 1) ||
				 (dxi == 1 && position[dxi] == 0));

		position[dxi] += 1 - 2 * dxi;
		//fprintf( stderr, "(%i, %i, %i)\n", position[0], position[1], position[2]);

		// new vessel node
		last = next;
		next = vesselGraph->octree->at(position[0], position[1], position[2]);
		if (next == 0) {
			next = new VesselNode(position[0], position[1], position[2]);
			next->setType(VESSEL);
			vesselGraph->addVesselNode(next);
		}
		vesselGraph->addVesselSegment(new VesselSegment(last, next));
		vesselGraph->vesselSegments[vesselGraph->countVesselSegments - 1]->radius
				= 25;

		position[dxi] += 1 - 2 * dxi;
		last = next;
		next = vesselGraph->octree->at(position[0], position[1], position[2]);
		if (next == 0) {
			next = new VesselNode(position[0], position[1], position[2]);
			next->setType(VESSEL);
			vesselGraph->addVesselNode(next);
		}
		vesselGraph->addVesselSegment(new VesselSegment(last, next));
		vesselGraph->vesselSegments[vesselGraph->countVesselSegments - 1]->radius
				= 25;
		i++;

		//fprintf( stdout, "%i %i %i %i %i %i\n",
		//		last->position[0], last->position[1], last->position[2],
		//		next->position[0]-last->position[0], next->position[1]-last->position[1], next->position[2]-last->position[2]);
	}
	next->pressure = 0;
	next->setType(ROOT);
}

void updateOxygen(float **oxygen, float **oxygen2, VesselGraph *vesselGraph,
		int maxit, int MaxError) {
	float **temp;
	float maxError = 0;
	float lastmaxError = 0;
	int it = 0;
	do {
		lastmaxError = maxError;
		maxError = 0;
		//for( int it=0; it<1000; it++){

		//for( int i=0; i<vesselGraph->countVesselNodes; i++)
		//	oxygen[vesselGraph->vesselNodes[i]->position[0]][vesselGraph->vesselNodes[i]->position[1]] = 1.;

		// CALCULATE OXYGEN CONCENTRATION
		//float D = 1.;
		float dt = 10.;
		float dx = 10.;
		for (int x = 0; x < DOMAIN_SIZE_X; x++) {
			for (int y = 0; y < DOMAIN_SIZE_Y; y++) {
				//fprintf( stderr, "%i %i\n", x,y);
				float diffusion = 0.;
				float consumption = 1. / 10000.; // consumption rate normal tissue
				//- 0.004*(1.-oxygen[x][y]); // source rate vessels
				for (int xx = x - 1; xx <= x + 1; xx += 2)
					for (int yy = y - 1; yy <= y + 1; yy += 2)
						if (xx >= 0 && xx < DOMAIN_SIZE_X && yy >= 0 && yy
								< DOMAIN_SIZE_Y) {
							//fprintf( stderr, "%i %i %i %i\n", x,y,xx,yy);
							diffusion += oxygen[xx][yy] - oxygen[x][y];
						}
				diffusion /= (dx * dx);
				oxygen2[x][y] = oxygen[x][y] + (diffusion - consumption) * dt;
				if (maxError < fabs(oxygen2[x][y] - oxygen[x][y]))
					maxError = fabs(oxygen2[x][y] - oxygen[x][y]);
			}
		}

		for (int i = 0; i < vesselGraph->countVesselNodes; i++)
			if (vesselGraph->vesselNodes[i]->getType() == VESSEL
					&& vesselGraph->vesselNodes[i]->countNeighbors > 0
					&& fabs(
							vesselGraph->vesselNodes[i]->pressure
									- vesselGraph->vesselNodes[i]->neighbors[0]->pressure)
							> 1e-4)
				oxygen2[(int) vesselGraph->vesselNodes[i]->position[0]][(int) vesselGraph->vesselNodes[i]->position[1]]
						+= 0.004
								* (1.
										- oxygen[(int) vesselGraph->vesselNodes[i]->position[0]][(int) vesselGraph->vesselNodes[i]->position[1]])
								* dt; // source rate vessels ;

		for (int x = 0; x < DOMAIN_SIZE_X; x++)
			for (int y = 0; y < DOMAIN_SIZE_Y; y++)
				if (oxygen2[x][y] < 0.)
					oxygen2[x][y] = 0.;

		temp = oxygen;
		oxygen = oxygen2;
		oxygen2 = temp;

		//for( int i=0; i<vesselGraph->countVesselNodes; i++)
		//	oxygen[vesselGraph->vesselNodes[i]->position[0]][vesselGraph->vesselNodes[i]->position[1]] = 1.;
		it++;
		fprintf(stderr, "\rmax error = %lf \b", maxError);
	} while ((maxError > MaxError && fabs(maxError - lastmaxError) > 1e-20)
			&& it < maxit);
	fprintf(stderr, "INFO: Oxygen Diffusion finished after %i iterations\n", it);

}

void updateVEGF(float **oxygen, float **vegf) {
	for (int x = 0; x < DOMAIN_SIZE_X; x++)
		for (int y = 0; y < DOMAIN_SIZE_Y; y++)
			vegf[x][y] = 0.;

	int dx = 10;
	int Rvegf = 200 / dx;
	float fRvegf = (float) Rvegf;
	//int dist = Rvegf / dx;
	for (int x = 0; x < DOMAIN_SIZE_X; x++)
		for (int y = 0; y < DOMAIN_SIZE_Y; y++)
			if (oxygen[x][y] < 0.1)
				for (int xx = x - Rvegf; xx <= x + Rvegf; xx++)
					for (int yy = y - Rvegf; yy <= y + Rvegf; yy++)
						if (xx >= 0 && xx < DOMAIN_SIZE_X && yy >= 0 && yy
								< DOMAIN_SIZE_Y) {
							float fdist = sqrt(pow(xx - x, 2) + pow(yy - y, 2));
							if (fdist < fRvegf)
								vegf[xx][yy] += 1 - fdist / fRvegf;
						}

}

void setInterTipConnections(VesselGraph *vesselGraph,
		int &countInterTipConnections) {
	//int countPreviousVesselNodes = vesselGraph->countVesselNodes;
	countInterTipConnections = 0;
	for (int v = 0; v < vesselGraph->countVesselNodes; v++)
		if (vesselGraph->vesselNodes[v]->getType() == TIP //|| vesselGraph->vesselNodes[v]->getType() == ROOT
				) {
			float x = vesselGraph->vesselNodes[v]->position[0];
			float y = vesselGraph->vesselNodes[v]->position[1];
			float z = vesselGraph->vesselNodes[v]->position[2];

			VesselNode *vn;
			//fprintf( stderr, "(%.0f,%.0f,%.0f)\n",x,y,z);
			//if(x >= BRANCH_LENGTH && vesselGraph->octree->at(x - BRANCH_LENGTH, y, z)!=0)
				//fprintf( stderr, "->(%.0f,%.0f,%.0f)\n",vesselGraph->octree->at(x - BRANCH_LENGTH, y, z)->position[0],y,z);

			if (x >= BRANCH_LENGTH) {
				vn = vesselGraph->octree->at(x - BRANCH_LENGTH, y, z);
				if (vn != 0 && (vn->getType() == TIP //|| vn->getType() == ROOT
						)) {
					//fprintf( stderr, "Add Inter Tip Connection (x-direction) between (%.0f,%.0f,%.0f) and (%.0f,%.0f,%.0f)\n",
					//		x,y,z, x - BRANCH_LENGTH, y, z);
					countInterTipConnections++;
					vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
							vesselGraph->vesselNodes[v], vn));
				}
			}

			if (y >= BRANCH_LENGTH) {
				vn = vesselGraph->octree->at(x, y - BRANCH_LENGTH, z);
				if (vn != 0 && (vn->getType() == TIP //|| vn->getType() == ROOT
						)) {
					//fprintf( stderr, "Add Inter Tip Connection\n");
					countInterTipConnections++;
					vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
							vesselGraph->vesselNodes[v], vn));
				}
			}

			if (z >= BRANCH_LENGTH) {
				vn = vesselGraph->octree->at(x, y, z - BRANCH_LENGTH);
				if (vn != 0 && (vn->getType() == TIP //|| vn->getType() == ROOT
						)) {
					//fprintf( stderr, "Add Inter Tip Connection\n");
					countInterTipConnections++;
					vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
							vesselGraph->vesselNodes[v], vn));
				}
			}

		}
}

void setInterTipConnectionsAll(VesselGraph *vesselGraph,
		int &countInterTipConnections) {
	//int countPreviousVesselNodes = vesselGraph->countVesselNodes;
	countInterTipConnections = 0;
	for (int v = 0; v < vesselGraph->countVesselNodes; v++)
		if (vesselGraph->vesselNodes[v]->getType() == TIP) {
			float x = vesselGraph->vesselNodes[v]->position[0];
			float y = vesselGraph->vesselNodes[v]->position[1];
			float z = vesselGraph->vesselNodes[v]->position[2];

			VesselNode *vn;

			if (x >= BRANCH_LENGTH) {
				vn = vesselGraph->octree->at(x - BRANCH_LENGTH, y, z);
				if (vn != 0 && vn->getType() != ROOT) {
					//fprintf( stderr, "Add Inter Tip Connection\n");
					countInterTipConnections++;
					vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
							vesselGraph->vesselNodes[v], vn));
				}
			}

			if (y >= BRANCH_LENGTH) {
				vn = vesselGraph->octree->at(x, y - BRANCH_LENGTH, z);
				if (vn != 0 && vn->getType() != ROOT) {
					//fprintf( stderr, "Add Inter Tip Connection\n");
					countInterTipConnections++;
					vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
							vesselGraph->vesselNodes[v], vn));
				}
			}

			if (z >= BRANCH_LENGTH) {
				vn = vesselGraph->octree->at(x, y, z - BRANCH_LENGTH);
				if (vn != 0 && vn->getType() != ROOT) {
					//fprintf( stderr, "Add Inter Tip Connection\n");
					countInterTipConnections++;
					vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment(
							vesselGraph->vesselNodes[v], vn));
				}
			}

		}
}

void removeInterTipConnections(VesselGraph *vesselGraph,
		int &countInterTipConnections) {
	//int countPreviousVesselNodes = vesselGraph->countVesselNodes;
	countInterTipConnections = 0;
	for (int v = vesselGraph->countVesselSegments - 1; v >= 0; v--)
		if( vesselGraph->vesselSegments[v]->vesselNodes[0]->getType() == TIP &&
			vesselGraph->vesselSegments[v]->vesselNodes[1]->getType() == TIP) {
			vesselGraph->removeVesselSegmentAndNeighbors(vesselGraph->vesselSegments[v]);
			delete vesselGraph->vesselSegments[v];
			countInterTipConnections++;
		}
}

void updateVesselGraph(float **oxygen, float **vegf, VesselGraph *vesselGraph,
		float dt, float time) {
	int countPreviousVesselNodes = vesselGraph->countVesselNodes;
	int countPreviousVesselSegments = vesselGraph->countVesselSegments;

	for (int i = 0; i < countPreviousVesselNodes; i++) {

		if (time - vesselGraph->vesselNodes[i]->time >= 100)
			switch (vesselGraph->vesselNodes[i]->getType()) {
			case SPROUT:
			case TIP:
				vesselGraph->vesselNodes[i]->setType(VESSEL);
				break;
			}

		switch (vesselGraph->vesselNodes[i]->getType()) {
		case ROOT: //break;

		case VESSEL:
			// sprout initiation
		case SPROUT:
			// sprout migration
			if ( (vesselGraph->vesselNodes[i]->countNeighbors == 0
				  || (vesselGraph->vesselNodes[i]->countNeighbors < 3
							&& vesselGraph->vesselNodes[i]->neighbors[0]->countNeighbors < 3
							&& (vesselGraph->vesselNodes[i]->countNeighbors	== 1 || vesselGraph->vesselNodes[i]->neighbors[1]->countNeighbors < 3)
							&& vegf[(int) vesselGraph->vesselNodes[i]->position[0]][(int) vesselGraph->vesselNodes[i]->position[1]] > 0.))
							&& dt / 5. >= RAND01) {
				int candidateX = -1;
				int candidateY = -1;
				int x = vesselGraph->vesselNodes[i]->position[0];
				int y = vesselGraph->vesselNodes[i]->position[1];
				float maxVEGF = vegf[x][y];

				for (int xx = x - 1; xx <= x + 1; xx += 2)
					if (xx >= 0 && xx < DOMAIN_SIZE_X /*&& yy>=0 && yy<DOMAIN_SIZE*/
					&& vegf[xx][y] > maxVEGF && (y == 0
							|| vesselGraph->octree->at(xx, y - 1, 0) == 0)
							&& (vesselGraph->octree->at(xx, y, 0) == 0)
							&& (vesselGraph->octree->at(xx, y + 1, 0) == 0)) {
						maxVEGF = vegf[xx][y];
						candidateX = xx;
						candidateY = y;
					}

				for (int yy = y - 1; yy <= y + 1; yy += 2)
					if ( /*xx>=0 && xx<DOMAIN_SIZE &&*/yy >= 0 && yy
							< DOMAIN_SIZE_Y && vegf[x][yy] > maxVEGF && (x == 0
							|| vesselGraph->octree->at(x - 1, yy, 0) == 0)
							&& vesselGraph->octree->at(x, yy, 0) == 0
							&& vesselGraph->octree->at(x + 1, yy, 0) == 0) {
						maxVEGF = vegf[x][yy];
						candidateX = x;
						candidateY = yy;
					}
				//fprintf( stderr, "Sprout! maxVEGF=%e\n", maxVEGF);

				if (candidateX != -1) {
					VesselNode *next =
							new VesselNode(candidateX, candidateY, 0);
					next->setType(TIP);
					vesselGraph->addVesselNode(next);
					vesselGraph->addVesselSegment(new VesselSegment(
							vesselGraph->vesselNodes[i], next));
					if (vesselGraph->vesselNodes[i]->getType() == SPROUT)
						next->time = vesselGraph->vesselNodes[i]->time;
					else
						next->time = time;
				}

			}
			break;

		case TIP:

			// new position
			int pos[3];
			for (int d = 0; d < 2; d++) {
				pos[d]
						= 2 * vesselGraph->vesselNodes[i]->position[d]
								- vesselGraph->vesselNodes[i]->neighbors[0]->position[d];
			}

			if (pos[0] >= 0 && pos[0] < DOMAIN_SIZE_X && pos[1] >= 0 && pos[1]
					< DOMAIN_SIZE_Z && dt / 5. >= RAND01) {
				VesselNode *next = vesselGraph->octree->at(pos[0], pos[1], 0);
				vesselGraph->vesselNodes[i]->setType(SPROUT);
				if (next == 0) {
					//fprintf( stderr, "Sprout Migr.!\n");

					next = new VesselNode(pos[0], pos[1], 0);
					next->setType(TIP);
					next->time = vesselGraph->vesselNodes[i]->time;
					vesselGraph->addVesselNode(next);
					//if( time - vesselGraph->vesselNodes[i]->time >= 100)
					//next->setType( VESSEL);
				}//else fprintf( stderr, "Sprout meets Network !\n");
				vesselGraph->addVesselSegment(new VesselSegment(
						vesselGraph->vesselNodes[i], next));
			}

			// termination
			/*if( time - vesselGraph->vesselNodes[i]->time >= 100)
			 {
			 VesselNode *next = vesselGraph->vesselNodes[i];
			 while( next->time == vesselGraph->vesselNodes[i]->time){
			 next->setType( VESSEL);
			 if( next->neighbors[0]->type == SPROUT && next->neighbors[0]->time == vesselGraph->vesselNodes[i]->time)
			 next = next->neighbors[0];
			 else
			 next = next->neighbors[1];
			 }
			 }*/

			break;

		default:
			break;
		}
		//if( vegf[vesselGraph->vesselNodes[i]->position[0]][vesselGraph->vesselNodes[i]->position[1]] > 0.
		//		&& 1./50. >= RAND01){
	}

	for (int i = 0; i < countPreviousVesselSegments; i++) {
		if ((vesselGraph->vesselSegments[i]->vesselNodes[0]->getType()
				== VESSEL
				|| vesselGraph->vesselSegments[i]->vesselNodes[0]->getType()
						== ROOT)
				&& (vesselGraph->vesselSegments[i]->vesselNodes[1]->getType()
						== VESSEL
						|| vesselGraph->vesselSegments[i]->vesselNodes[1]->getType()
								== ROOT)) {
			if (time - vesselGraph->vesselSegments[i]->vesselNodes[0]->time
					>= 124
					&& vegf[(int) vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]][(int) vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
							> 0 && vesselGraph->vesselSegments[i]->radius < 25)
				vesselGraph->vesselSegments[i]->radius += 0.4 * dt;

			if (vesselGraph->vesselSegments[i]->shear < 0.002) {
				if (dt / 20. >= RAND01 /*&& vesselGraph->vesselSegments[i]->radius < 20*/) {
					fprintf(stderr, "Vessel Regression!\n");
					VesselSegment *removedVS = vesselGraph->vesselSegments[i];
					VesselNode *removedVN0 = removedVS->vesselNodes[0];
					VesselNode *removedVN1 = removedVS->vesselNodes[1];

					//removedVS->vesselNodes[0]->removeNeighbor( removedVS->vesselNodes[1]);
					//removedVS->vesselNodes[1]->removeNeighbor( removedVS->vesselNodes[0]);

					vesselGraph->removeVesselSegment(removedVS);
					delete removedVS;

					if (removedVN0->countNeighbors == 0
							&& removedVN0->getType() != ROOT) {
						VesselNode *removedVN = removedVN0;
						vesselGraph->removeVesselNode(removedVN);
						delete removedVN;
					}
					if (removedVN1->countNeighbors == 0
							&& removedVN1->getType() != ROOT) {
						VesselNode *removedVN = removedVN1;
						vesselGraph->removeVesselNode(removedVN);
						delete removedVN;
					}
				}
			}
			//else

		}
	}
}

int perfusion(int argc, char** argv, VesselGraph *vesselGraph) {



	//bool vessels[DOMAIN_SIZE][DOMAIN_SIZE][DOMAIN_SIZE];
	if(vesselGraph == 0)
		vesselGraph = new VesselGraph();
	float time = 0.;
	float minShear=0, maxShear=0;
	//vesselGraph->addVesselSegment(
	//		new VesselSegment( new VesselNode(0,0,0), new VesselNode(1,0,0)));
	printColorMapLegend( (char*)"legend.eps", (char*)"1.1m", (char*)"12.12m");



	fprintf( stderr, "\n\nINITIALIZE BLOOD VESSEL NETWORK\n\n");

	setInitialVesselGraph(vesselGraph);
	//vesselGraph->updatePressure(0,0,0);
	vesselGraph->updateRadius();


	//vesselGraph->printToPovray( "before.pov", 0);

	/*for( double t=0.; t<20.; t=t+0.1){
		printf( "%lf %lf\n", t, VesselGraph::AIF::ArterialInputFunction( t, 10, VesselGraph::AIF::PARKER));
	}
	exit(0);*/
	{
	int countInterTipConnections = 0;
	setInterTipConnections(vesselGraph, countInterTipConnections);
	vesselGraph->updatePressureNEW2();
	vesselGraph->updateFlow();
	vesselGraph->updateShear();
	fprintf(stderr, "countInterTipConnections=%i\n",countInterTipConnections);
	//vesselGraph->printToEPS( "beforeRemodelling.eps", 0, NULL);
	removeInterTipConnections(vesselGraph, countInterTipConnections);
	}
//exit(0);

	fprintf( stderr, "\n\nPERFORM SHEAR STRESS HOMOGENIZATION\n\n");


	for (int it = 0; it < 100; it++) {
		int countInterTipConnections = 0;
		vesselGraph->updateSegmentTypes();
		setInterTipConnections(vesselGraph, countInterTipConnections);
		fprintf(stderr, "countInterTipConnections=%i\n",
				countInterTipConnections);
		//vesselGraph->printToPovray((char*) "before.pov", 0);
		vesselGraph->updateRadius();
		//vesselGraph->updatePressure( 0,0,0);
		vesselGraph->updatePressureNEW2();
		vesselGraph->updateFlow();
		vesselGraph->updateShear();
		/*if(it==0){
		vesselGraph->printToPovray((char*) "before.pov", 0);
		FILE *fp = fopen("vesselBefore.dat", "w+");
		for (int i = 0; i < vesselGraph->countVesselSegments; i++)
			for (int ii = 0; ii < 2; ii++)
				fprintf(
						fp,
						"%i %i %e %e %e %i %lf %lf %lf\n",
						(int) vesselGraph->vesselSegments[i]->vesselNodes[ii]->position[0],
						(int) vesselGraph->vesselSegments[i]->vesselNodes[ii]->position[1],
						vesselGraph->vesselSegments[i]->flow,
						vesselGraph->vesselSegments[i]->shear,
						vesselGraph->vesselSegments[i]->vesselNodes[ii]->pressure,
						vesselGraph->vesselSegments[i]->vesselNodes[ii]->getType(),
						vesselGraph->vesselSegments[i]->radius,
						(vesselGraph->vesselSegments[i]->shear - minShear)
								/ (maxShear - minShear), (maxShear
								- vesselGraph->vesselSegments[i]->shear)
								/ (maxShear - minShear));
		fclose(fp);

		}*/

		/*char filename[512];
		sprintf( filename, "remodelling%i.eps", it);
		vesselGraph->printToEPS( filename, 0, NULL);
		 */
		/*sprintf( filename, "remodelling%i.pov", it);
		vesselGraph->printToPovray( filename, 0);*/

		// STATISTICS
		vesselGraph->writeStatisticsToFile( "statisticsPressure.dat", 0);
		vesselGraph->writeStatisticsToFile( "statisticsFlow.dat", 1);
		vesselGraph->writeStatisticsToFile( "statisticsVelocity.dat", 2);
		vesselGraph->writeStatisticsToFile( "statisticsShearStress.dat", 3);
		// END STATISTICS


		removeInterTipConnections(vesselGraph, countInterTipConnections);
		fprintf(stderr, "\ncountInterTipConnections=%i removed!\n",
				countInterTipConnections);

		// REMODELING
		// [min,max] shear stress

		minShear = maxShear = vesselGraph->vesselSegments[0]->shear;
		for (int i = 1; i < vesselGraph->countVesselSegments; i++)
			//if (vesselGraph->vesselSegments[i]->vesselNodes[0]->getType() == TIP
			//		|| vesselGraph->vesselSegments[i]->vesselNodes[1]->getType() == TIP)
			{
				if (minShear > vesselGraph->vesselSegments[i]->shear)
					minShear = vesselGraph->vesselSegments[i]->shear;
				if (maxShear < vesselGraph->vesselSegments[i]->shear)
					maxShear = vesselGraph->vesselSegments[i]->shear;
			}

		// DEGENERATION
		float randShear;// = RAND01 * (maxShear - minShear) + minShear;
		for( int t=0; t<2; t++){
		for (int i = vesselGraph->countVesselSegments - 1; i >= 0; i--){
			randShear = RAND01 * (maxShear - minShear) + minShear;
			//if(1./4. > RAND01)
			if((vesselGraph->vesselSegments[i]->vesselNodes[0]->getType() == TIP &&
				vesselGraph->vesselSegments[i]->vesselNodes[0]->countNeighbors == 1) ||
			   (vesselGraph->vesselSegments[i]->vesselNodes[1]->getType() == TIP &&
				vesselGraph->vesselSegments[i]->vesselNodes[1]->countNeighbors == 1)) {

				if (vesselGraph->vesselSegments[i]->index != i) {
					fprintf(stderr, "WORNG INDEX!\n");
					exit(0);
				}
				if (vesselGraph->vesselSegments[i]->shear < randShear) {
					if (vesselGraph->vesselSegments[i]->vesselNodes[0]->getType() != ROOT && vesselGraph->vesselSegments[i]->vesselNodes[0]->countNeighbors==2)
						vesselGraph->vesselSegments[i]->vesselNodes[0]->setType(TIP);
					if (vesselGraph->vesselSegments[i]->vesselNodes[1]->getType()!= ROOT && vesselGraph->vesselSegments[i]->vesselNodes[1]->countNeighbors==2)
						vesselGraph->vesselSegments[i]->vesselNodes[1]->setType(TIP);
					VesselNode *rvn0 = vesselGraph->vesselSegments[i]->vesselNodes[0],
							   *rvn1 = vesselGraph->vesselSegments[i]->vesselNodes[1];
					vesselGraph->removeVesselSegmentAndNeighbors(vesselGraph->vesselSegments[i]);
					if( rvn0->countNeighbors==0)
						vesselGraph->removeVesselNode( rvn0);
					if( rvn1->countNeighbors==0)
						vesselGraph->removeVesselNode( rvn1);
				}
			}
		}
		}

		// SPROUTING
		for( int t=0; t<10; t++){
		int countOldVesselNodes = vesselGraph->countVesselNodes;
		for (int i = 0; i < countOldVesselNodes; i++){
			randShear = RAND01 * (maxShear - minShear) + minShear;
			if( vesselGraph->vesselNodes[i]->countNeighbors <= 2 && vesselGraph->vesselNodes[i]->getType()!=ROOT)
			//if(vesselGraph->vesselSegments[i]->shear > randShear)
			{
				float sproutingProb = 1/2.;//1./2.;//1./5.
				float dx[6][3] = {{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0}};
				//float dx[6][3] = {0,0,0};
				int count = 0;
				if( (int)vesselGraph->vesselNodes[i]->position[0] < DOMAIN_SIZE_X-1 &&
						vesselGraph->octree->at(
						(int)vesselGraph->vesselNodes[i]->position[0]+1,
						(int)vesselGraph->vesselNodes[i]->position[1],
						(int)vesselGraph->vesselNodes[i]->position[2]) == 0 && RAND01 < sproutingProb
						){
					//fprintf(stderr, "CANDIDATE\n");
					dx[count][0]=1; dx[count][1]=0; dx[count][2]=0;
					count++;
					//dx[0]=1; dx[1]=0; dx[2]=0;
				}randShear = RAND01 * (maxShear - minShear) + minShear;
				if( (int)vesselGraph->vesselNodes[i]->position[0] > 0 &&
						vesselGraph->octree->at(
								(int)vesselGraph->vesselNodes[i]->position[0]-1,
								(int)vesselGraph->vesselNodes[i]->position[1],
								(int)vesselGraph->vesselNodes[i]->position[2]) == 0 && RAND01 < sproutingProb
								){
					//fprintf(stderr, "CANDIDATE\n");
					dx[count][0]=-1; dx[count][1]=0; dx[count][2]=0;
					count++;
					//dx[0]=-1; dx[1]=0; dx[2]=0;
				}randShear = RAND01 * (maxShear - minShear) + minShear;
				if( (int)vesselGraph->vesselNodes[i]->position[1] < DOMAIN_SIZE_Y-1 &&
						vesselGraph->octree->at(
								(int)vesselGraph->vesselNodes[i]->position[0],
								(int)vesselGraph->vesselNodes[i]->position[1]+1,
								(int)vesselGraph->vesselNodes[i]->position[2]) == 0 && RAND01 < sproutingProb
								){
					//fprintf(stderr, "CANDIDATE\n");
					dx[count][0]=0; dx[count][1]=1; dx[count][2]=0;
					count++;
					//dx[0]=0; dx[1]=1; dx[2]=0;
				}randShear = RAND01 * (maxShear - minShear) + minShear;
				if( (int)vesselGraph->vesselNodes[i]->position[1] > 0 &&
						vesselGraph->octree->at(
								(int)vesselGraph->vesselNodes[i]->position[0],
								(int)vesselGraph->vesselNodes[i]->position[1]-1,
								(int)vesselGraph->vesselNodes[i]->position[2]) == 0 && RAND01 < sproutingProb
								){
					//fprintf(stderr, "CANDIDATE\n");
					dx[count][0]=0; dx[count][1]=-1; dx[count][2]=0;
					count++;
					//dx[0]=0; dx[1]=-1; dx[2]=0;
				}randShear = RAND01 * (maxShear - minShear) + minShear;
				if( (int)vesselGraph->vesselNodes[i]->position[2] < DOMAIN_SIZE_Z-1 &&
						vesselGraph->octree->at(
								(int)vesselGraph->vesselNodes[i]->position[0],
								(int)vesselGraph->vesselNodes[i]->position[1],
								(int)vesselGraph->vesselNodes[i]->position[2]+1) == 0 && RAND01 < sproutingProb
								){
					//fprintf(stderr, "CANDIDATE\n");
					dx[count][0]=0; dx[count][1]=0; dx[count][2]=1;
					count++;
					//dx[0]=0; dx[1]=0; dx[2]=1;
				}
				if( (int)vesselGraph->vesselNodes[i]->position[2] > 0 &&
						vesselGraph->octree->at(
								(int)vesselGraph->vesselNodes[i]->position[0],
								(int)vesselGraph->vesselNodes[i]->position[1],
								(int)vesselGraph->vesselNodes[i]->position[2]-1) == 0 && RAND01 < sproutingProb
								){
					//fprintf(stderr, "CANDIDATE\n");
					dx[count][0]=0; dx[count][1]=0; dx[count][2]=-1;
					count++;
					//dx[0]=0; dx[1]=0; dx[2]=-1;
				}

				int which = (int) (RAND01 * (float)count);
				if( dx[which][0]!=0 || dx[which][1]!=0 || dx[which][2]!=0){
					//fprintf(stderr, "SPROUT!\n");
					VesselNode *next = new VesselNode(vesselGraph->vesselNodes[i]->position[0]+dx[which][0], vesselGraph->vesselNodes[i]->position[1]+dx[which][1], vesselGraph->vesselNodes[i]->position[2]+dx[which][2]);
					if( vesselGraph->vesselNodes[i]->getType()!=ROOT )
						vesselGraph->vesselNodes[i]->setType( VESSEL);
					next->setType(TIP);
					vesselGraph->addVesselNode(next);
					vesselGraph->addVesselSegmentAndNeighbors(new VesselSegment( vesselGraph->vesselNodes[i], next));

				}
			}
		}
		}
		/*for (int i = 0; i < vesselGraph->countVesselNodes; i++){
			if( vesselGraph->vesselNodes[i]->getType()==TIP && vesselGraph->vesselNodes[i]->countNeighbors!=1){
				fprintf(stderr, "TIP ERROR!\n");
				exit(0);
			}
			if( vesselGraph->vesselNodes[i]->countNeighbors>3){
				fprintf(stderr, "Toomany neighbors!\n");
				exit(0);
			}
		}*/
	}

	{
		int countInterTipConnections = 0;
		setInterTipConnections(vesselGraph, countInterTipConnections);
		fprintf(stderr, "countInterTipConnections=%i\n",
				countInterTipConnections);

		vesselGraph->updateRadius();
		vesselGraph->updatePressureNEW();

		char filename[512] = "afterRemodelling.pov";
		vesselGraph->printToPovray( filename, 0);
		vesselGraph->printToEPS( "afterRemodelling.eps", 0, NULL);
	}
	//exit( 0);

	/*float noise = 0.5;
	 for( int i=0; i<vesselGraph->countVesselNodes; i++){
	 vesselGraph->vesselNodes[i]->position[0] += noise *(0.5+RAND01);
	 vesselGraph->vesselNodes[i]->position[1] += noise *(0.5+RAND01);
	 vesselGraph->vesselNodes[i]->position[2] += noise *(0.5+RAND01);
	 }*/

	vesselGraph->updateRadius();
	vesselGraph->updateFlow();
	//for(int vs=0; vs<vesselGraph->countVesselSegments; vs++)
	//	vesselGraph->vesselSegments[vs]->flow *= 100000. / fabs(vesselGraph->vesselSegments[vs]->flow);
	vesselGraph->updateShear();

	// SET PERMEABILITY
	for( int i=0; i<vesselGraph->countVesselSegments; i++){
		// Tumor in (x-20)^2 + (y-20)^2 <= 10^2
		float x = 0.5*(vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0] + vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]);
		float y = 0.5*(vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1] + vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]);
		if( pow(x-25,2) + pow(y-25,2) <= pow(20,2))
			// leaking vessels in tumor
			vesselGraph->vesselSegments[i]->permeability = 100;
		else
			// functional vascularisation in surrounding tissue
			vesselGraph->vesselSegments[i]->permeability = 10;
	}

	/*int countInterTipConnections = 0;
	setInterTipConnections(vesselGraph, countInterTipConnections);
	vesselGraph->updatePressureNEW();
	removeInterTipConnections(vesselGraph, countInterTipConnections);
	vesselGraph->updateShear();*/
	//vesselGraph->updateFlow();
	//vesselGraph->printToPovray((char*) "test.pov", 0);
	//vesselGraph->printToPovray((char*) "test.pov", 0);

	/*FILE *fp = fopen("vessel.dat", "w+");
	for (int i = 0; i < vesselGraph->countVesselSegments; i++)
		for (int ii = 0; ii < 2; ii++)
			fprintf(
					fp,
					"%i %i %e %e %e %i %lf %lf %lf\n",
					(int) vesselGraph->vesselSegments[i]->vesselNodes[ii]->position[0],
					(int) vesselGraph->vesselSegments[i]->vesselNodes[ii]->position[1],
					vesselGraph->vesselSegments[i]->flow,
					vesselGraph->vesselSegments[i]->shear,
					vesselGraph->vesselSegments[i]->vesselNodes[ii]->pressure,
					vesselGraph->vesselSegments[i]->vesselNodes[ii]->getType(),
					vesselGraph->vesselSegments[i]->radius,
					(vesselGraph->vesselSegments[i]->shear - minShear)
							/ (maxShear - minShear), (maxShear
							- vesselGraph->vesselSegments[i]->shear)
							/ (maxShear - minShear));
	fclose(fp);*/

	// OUTPUT EXACT

	double mean_volume_fraction_vessels = 0.;
	double mean_flow = 0.;
	double mean_flow_per_vessel_volume = 0.;
	double mean_surface = 0.;
	double mean_surface_per_vessel_volume = 0.;

	FILE *fp = fopen("vessel.dat", "w+");
	for (int i = 0; i < vesselGraph->countVesselNodes; i++){
		// Average Flow
		double flow = 0.;
		for( int v=0; v<vesselGraph->vesselNodes[i]->countNeighbors; v++)
		{
			// out flow
			if( vesselGraph->vesselNodes[i]->pressure >= vesselGraph->vesselNodes[i]->neighbors[v]->pressure)
				flow += fabs(vesselGraph->vesselNodes[i]->branches[v]->flow);
		}

		// Average Shear Stress
		double shear = 0.;
		for( int v=0; v<vesselGraph->vesselNodes[i]->countNeighbors; v++)
		{
			// out flow
			//if( vesselGraph->vesselNodes[i]->pressure >= vesselGraph->vesselNodes[i]->neighbors[v]->pressure)
				shear += fabs(vesselGraph->vesselNodes[i]->branches[v]->shear);
		}

		// Volume Fraction (Vessels)
		double volume_total = LATTICE_CONSTANT*LATTICE_CONSTANT*LATTICE_CONSTANT;
		double volume_vessels = 0.;
		for( int v=0; v<vesselGraph->vesselNodes[i]->countNeighbors; v++)
			volume_vessels += 0.5*LATTICE_CONSTANT
			               * PI * vesselGraph->vesselNodes[i]->branches[v]->radius*vesselGraph->vesselNodes[i]->branches[v]->radius;

		// Surface of Vessels
		double surface_vessels = 0.;
		for( int v=0; v<vesselGraph->vesselNodes[i]->countNeighbors; v++)
			surface_vessels += LATTICE_CONSTANT * PI * vesselGraph->vesselNodes[i]->branches[v]->radius;


		fprintf(
					fp,
					"%i %i %i  %e %e %e %e %e\n",

					// Voxel
					(int) vesselGraph->vesselNodes[i]->position[0],
					(int) vesselGraph->vesselNodes[i]->position[1],
					(int) vesselGraph->vesselNodes[i]->position[2],

					// Volume Fractions
					volume_vessels/volume_total,

					// Vessel Surfaces
					surface_vessels,

					// Flow
					flow,

					// Shear Stress
					shear,

					// Pressure
					vesselGraph->vesselNodes[i]->pressure);

		mean_volume_fraction_vessels += volume_vessels / volume_total;
		mean_flow += flow;
		mean_flow_per_vessel_volume += flow / volume_vessels;
		mean_surface = surface_vessels;
		mean_surface_per_vessel_volume = surface_vessels / volume_vessels;

	}
	fclose(fp);

	double volume_total_domain = DOMAIN_SIZE_X*DOMAIN_SIZE_Y*DOMAIN_SIZE_Z;//*LATTICE_CONSTANT*LATTICE_CONSTANT*LATTICE_CONSTANT;
	fprintf( stderr, "mean_volume_fraction_vessels = %lf\n", mean_volume_fraction_vessels/volume_total_domain);
	fprintf( stderr, "mean_flow = %lf\n", mean_flow/volume_total_domain);
	fprintf( stderr, "mean_flow_per_vessel_volume = %lf\n", mean_flow_per_vessel_volume/volume_total_domain);
	fprintf( stderr, "mean_surface = %lf\n", mean_surface/volume_total_domain);
	fprintf( stderr, "mean_surface_per_vessel_volume = %lf\n", mean_surface_per_vessel_volume/volume_total_domain);

	printFlowMap( (char*)"flow.60x60x60.eps", vesselGraph, 60.,60.,60.);
	printBloodVolumeFractionMap( (char*)"f_p.60x60x60.eps", vesselGraph, 60.,60.,60.);
	printInterstitialSpaceVolumeFractionMap( (char*)"f_i.60x60x60.eps", vesselGraph, 60.,60.,60.);
	printPermeabilityMap( (char*)"Kps.60x60x60.eps", vesselGraph, 60.,60.,60.);

	printFlowMap( (char*)"flow.300x300x3000.eps", vesselGraph, 300.,300.,3000.);
	printBloodVolumeFractionMap( (char*)"f_p.300x300x3000.eps", vesselGraph, 300.,300.,3000.);
	printInterstitialSpaceVolumeFractionMap( (char*)"f_i.300x300x3000.eps", vesselGraph, 300.,300.,3000.);
	printPermeabilityMap( (char*)"Kps.300x300x3000.eps", vesselGraph, 300.,300.,3000.);
	//exit(0);
	// OUTPUT END

	// OUTPUT MAP

	fp = fopen("vessel.dat", "w+");
	for (int i = 0; i < vesselGraph->countVesselNodes; i++){
		// Average Flow
		double flow = 0.;
		for( int v=0; v<vesselGraph->vesselNodes[i]->countNeighbors; v++)
		{
			// out flow
			if( vesselGraph->vesselNodes[i]->pressure >= vesselGraph->vesselNodes[i]->neighbors[v]->pressure)
				flow += fabs(vesselGraph->vesselNodes[i]->branches[v]->flow);
		}

		// Average Shear Stress
		double shear = 0.;
		for( int v=0; v<vesselGraph->vesselNodes[i]->countNeighbors; v++)
		{
			// out flow
			//if( vesselGraph->vesselNodes[i]->pressure >= vesselGraph->vesselNodes[i]->neighbors[v]->pressure)
				shear += fabs(vesselGraph->vesselNodes[i]->branches[v]->shear);
		}

		// Volume Fraction (Vessels)
		double volume_total = LATTICE_CONSTANT*LATTICE_CONSTANT*LATTICE_CONSTANT;
		double volume_vessels = 0.;
		for( int v=0; v<vesselGraph->vesselNodes[i]->countNeighbors; v++)
			volume_vessels += 0.5*LATTICE_CONSTANT
			               * PI * vesselGraph->vesselNodes[i]->branches[v]->radius*vesselGraph->vesselNodes[i]->branches[v]->radius;

		// Surface of Vessels
		double surface_vessels = 0.;
		for( int v=0; v<vesselGraph->vesselNodes[i]->countNeighbors; v++)
			surface_vessels += LATTICE_CONSTANT * PI * vesselGraph->vesselNodes[i]->branches[v]->radius;


		fprintf(
					fp,
					"%i %i %i  %e %e %e %e %e\n",

					// Voxel
					(int) vesselGraph->vesselNodes[i]->position[0],
					(int) vesselGraph->vesselNodes[i]->position[1],
					(int) vesselGraph->vesselNodes[i]->position[2],

					// Volume Fractions
					volume_vessels/volume_total,

					// Vessel Surfaces
					surface_vessels,

					// Flow
					flow,

					// Shear Stress
					shear,

					// Pressure
					vesselGraph->vesselNodes[i]->pressure);

		mean_volume_fraction_vessels += volume_vessels / volume_total;
		mean_flow += flow;
		mean_flow_per_vessel_volume += flow / volume_vessels;
		mean_surface = surface_vessels;
		mean_surface_per_vessel_volume = surface_vessels / volume_vessels;

	}
	fclose(fp);

	//exit( 0);

	// OUTPUT END

	/*last = new VesselNode(100,0,0);
	 next = new VesselNode(99,0,0);
	 vesselGraph->addVesselNode( last);
	 vesselGraph->addVesselNode( next);
	 vesselGraph->addVesselSegment( new VesselSegment( last, next));
	 last->pressure = 0.; // kPa
	 last->setType( ROOT);
	 next->pressure = 0.; // kPa
	 next->setType( VESSEL);

	 last = new VesselNode(0,100,0);
	 next = new VesselNode(0,99,0);
	 vesselGraph->addVesselNode( last);
	 vesselGraph->addVesselNode( next);
	 vesselGraph->addVesselSegment( new VesselSegment( last, next));
	 last->pressure = 12.; // kPa
	 last->setType( ROOT);
	 next->pressure = 12.; // kPa
	 next->setType( VESSEL);*/



	// FLOW
	double **markerVessel = (double**) malloc(sizeof(double*) * DOMAIN_SIZE_X);
	double **markerVessel2 = (double**) malloc(sizeof(double*) * DOMAIN_SIZE_X);
	float ***markerTissue = (float***) malloc(sizeof(float**) * DOMAIN_SIZE_X);
	float ***markerTissue2 = (float***) malloc(sizeof(float**) * DOMAIN_SIZE_X);
	for (int x = 0; x < DOMAIN_SIZE_X; x++){
		markerVessel[x] = (double*) malloc(sizeof(double) * DOMAIN_SIZE_Y);
		markerVessel2[x] = (double*) malloc(sizeof(double) * DOMAIN_SIZE_Y);
		markerTissue[x] = (float**) malloc(sizeof(float*) * DOMAIN_SIZE_Y);
		markerTissue2[x] = (float**) malloc(sizeof(float*) * DOMAIN_SIZE_Y);
		for (int y = 0; y < DOMAIN_SIZE_Y; y++){
			markerVessel[x][y] = 0.;
			markerVessel2[x][y] = 0.;

			markerTissue[x][y] = (float*) malloc(sizeof(float) * DOMAIN_SIZE_Z);
			markerTissue2[x][y] = (float*) malloc(sizeof(float) * DOMAIN_SIZE_Z);
			for (int z = 0; z < DOMAIN_SIZE_Z; z++){
				markerTissue[x][y][z] = 0.;
				markerTissue2[x][y][z] = 0.;
			}
		}
	}
	double **before = markerVessel,
			**after = markerVessel2,
			**temp;
	/*int countInterTipConnections;
	setInterTipConnectionsAll(vesselGraph, countInterTipConnections);
	vesselGraph->updateRadius();
	vesselGraph->updatePressureNEW();
	//vesselGraph->updatePressure(0,0,0);
	vesselGraph->updateFlow();*/
	//vesselGraph->printToPovray( "after.pov", 0);


	/*vesselGraph->updateRadius();
	//vesselGraph->updatePressure( 0,0,0);
	vesselGraph->updatePressureNEW();
	vesselGraph->updateFlow();
	vesselGraph->updateShear();*/

	fprintf( stderr, "\n\nPERFORM CONTRAST AGENT PERFUSION\n\n");

	/*double time_step = 0.1; // sec
	double time_end  = 6000.; // sec
	double outputResolution = 10; //sec
	 */


	double time_step = 3.01592895 / 100.; // sec
	double time_end  = 1000.; // sec
	double outputResolution = 1.; //sec
	bool   outputAIF = true;

	float *b=0, *x=0, *dx=0;
	SparseMatrix<float> *A=0;

	fp = fopen("markerSpaTem.dat", "w");
	fclose(fp);
	double AIF = VesselGraph::AIF::ArterialInputFunction( 0., 0., VesselGraph::AIF::PARKER);
	{
		char filename[512];
		sprintf( filename, "markerMRI%05i.300x300x3000.eps", 0);
		printMarkerIntensityMap( filename, vesselGraph, markerTissue, 300,300,3000, outputAIF,VesselGraph::AIF::ArterialInputFunction( time / 60., 0., VesselGraph::AIF::PARKER));
		sprintf( filename, "markerMRI%05i.60x60x60.eps", 0);
		printMarkerIntensityMap( filename, vesselGraph, markerTissue, 60,60,60, outputAIF,VesselGraph::AIF::ArterialInputFunction( time / 60., 0., VesselGraph::AIF::PARKER));
		fprintf(stderr, "output: %i\n", (int)floor(time/outputResolution));
		fp = fopen("markerSpaTem.dat", "a+");
		for (int i = 0; i<vesselGraph->countVesselNodes; i++)
			fprintf(fp, "%f %f %e\n", vesselGraph->vesselNodes[i]->position[0], vesselGraph->vesselNodes[i]->position[1], vesselGraph->vesselNodes[i]->marker);
		fclose(fp);
	}
	for( double time=0.; time<=time_end; time+=time_step)
	{
		//vesselGraph->updateMarkerVesselsAndInterstitialSpace( 0.000001,0,0,markerTissue,markerTissue2);

		//vesselGraph->updateMarkerVesselsAndInterstitialSpaceNEW2( time_step,AIF,markerTissue,markerTissue2, A,b,x);
		vesselGraph->updateMarkerVesselsAndInterstitialSpaceNEW4( time_step,AIF,markerTissue,markerTissue2);

//->		vesselGraph->updateMarkerVesselsExplicit( time_step,AIF, x);
		//vesselGraph->updateMarkerVesselsExplicitAntiDiffusion( time_step,AIF, x);
		//vesselGraph->updateMarkerVesselsExplicit( time_step,AIF, dx,x);

		//vesselGraph->updateMarkerVesselsAndInterstitialSpaceBrix( dt,0,0,markerTissue,markerTissue2);
		//vesselGraph->updateMarkerVessels( 0.000001,0,0);
		/*fprintf(stderr, "%lf <> floor:%lf (%i) ceil:%lf\n",
				time/outputResolution,
				floor(time/outputResolution), time/outputResolution == floor(time/outputResolution),
				ceil(time/outputResolution));
		 */
		AIF = VesselGraph::AIF::ArterialInputFunction( (time+time_step) / 60., 0., VesselGraph::AIF::PARKER);

		if( (int)(time/outputResolution) != (int)((time+time_step)/outputResolution)){
			fprintf(stderr, "\rSimulation Time: %.0lfsec = %.1lfmin (%.1lf\%%)             \b", time, time/60, time/time_end*100.);
			//fprintf(stderr, "%lf <> floor:%lf ceil:%lf\n", time/outputResolution, floor(time/outputResolution), ceil(time/outputResolution));
			float ***temp = markerTissue; markerTissue = markerTissue2; markerTissue2 = temp;
			fp = fopen("marker.dat", "w+");
			for (int i = 0; i<vesselGraph->countVesselNodes; i++){
				fprintf(fp, "%f %f %e\n", vesselGraph->vesselNodes[i]->position[0], vesselGraph->vesselNodes[i]->position[1], vesselGraph->vesselNodes[i]->marker);
				//fprintf(stderr, "%e\n", vesselGraph->vesselNodes[i]->marker);
			}
			fclose(fp);
			char filename[512];
			/*sprintf( filename, "marker%i.pov", it/outputSkips);
			vesselGraph->printToPovray( filename, 1, markerTissue);
			//vesselGraph->printToPovray( filename, 1);
			 */
			/*sprintf( filename, "marker%i.eps", it/outputSkips);
			vesselGraph->printToEPS( filename, 1, markerTissue);
*/
			fprintf(stderr, "output: %i\n", (int)((time+time_step)/outputResolution));
			sprintf( filename, "markerMRI%05i.300x300x3000.eps", (int)((time+time_step)/outputResolution));
			printMarkerIntensityMap( filename, vesselGraph, markerTissue, 300,300,3000, outputAIF,AIF);
			sprintf( filename, "markerMRI%05i.60x60x60.eps", (int)((time+time_step)/outputResolution));
			printMarkerIntensityMap( filename, vesselGraph, markerTissue, 60,60,60, outputAIF,AIF);

			fp = fopen("markerSpaTem.dat", "a+");
			for (int i = 0; i<vesselGraph->countVesselNodes; i++)
				fprintf(fp, "%f %f %f %e %e\n", vesselGraph->vesselNodes[i]->position[0], vesselGraph->vesselNodes[i]->position[1], time, vesselGraph->vesselNodes[i]->marker,
						VesselGraph::AIF::ArterialInputFunction( (time/*+time_step*/
							- PI*pow(vesselGraph->vesselNodes[i]->branches[0]->radius,2)/fabs(vesselGraph->vesselNodes[i]->branches[0]->flow)*vesselGraph->vesselNodes[i]->position[0]*LATTICE_CONSTANT
						) / 60., 0., VesselGraph::AIF::PARKER));
			fclose(fp);
		}
	}
	free(x);
	free(b);
	delete A;

	//vesselGraph->printToPovray((char*) "before.pov");

	//removeInterTipConnections(vesselGraph, countInterTipConnections);
	//fprintf(stderr, "\ncountInterTipConnections=%i removed!\n",
	//		countInterTipConnections);

	// FREE MEMORY
	delete vesselGraph;

	for (int x = 0; x < DOMAIN_SIZE_X; x++){
		for (int y = 0; y < DOMAIN_SIZE_Y; y++){
			free( markerTissue[x][y]);
			free( markerTissue2[x][y]);
		}
		free( markerVessel[x]);
		free( markerVessel2[x]);
		free( markerTissue[x]);
		free( markerTissue2[x]);
	}
	free( markerVessel);
	free( markerVessel2);
	free( markerTissue);
	free( markerTissue2);


return 0;
exit( 0);

// INIT OXYGEN CONCENTRATION
float **oxygen = (float**) malloc(sizeof(float*) * DOMAIN_SIZE_X),
		**oxygen2 = (float**) malloc(sizeof(float*) * DOMAIN_SIZE_X);
for (int x = 0; x < DOMAIN_SIZE_X; x++) {
	oxygen[x] = (float*) malloc(sizeof(float) * DOMAIN_SIZE_Y);
	oxygen2[x] = (float*) malloc(sizeof(float) * DOMAIN_SIZE_Y);
	for (int y = 0; y < DOMAIN_SIZE_Y; y++)
		oxygen[x][y] = 0.;
}

updateOxygen(oxygen, oxygen2, vesselGraph, 10000, 1e-5);

float **vegf = (float**) malloc(sizeof(float*) * DOMAIN_SIZE_X);
for (int x = 0; x < DOMAIN_SIZE_X; x++)
	vegf[x] = (float*) malloc(sizeof(float) * DOMAIN_SIZE_Y);

updateVEGF(oxygen, vegf);

	for( int it=0; it<100000; it++){
		fprintf( stderr, "\r%i \b", it);
	for (int i = 0; i<vesselGraph->countVesselSegments; i++)
	if(vesselGraph->vesselSegments[i]->flow!=0){
		if( isnan(vesselGraph->vesselSegments[i]->flow)){
			fprintf(stderr, "SHIT!\n");
			exit( 0);
		}
		//fprintf(stderr, "YEAH! %e %e %e\n", before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
		//              				       [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
		//              			   * vesselGraph->vesselSegments[i]->flow * dt);

		/*if( vesselGraph->vesselSegments[i]->flow<0){
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
			     [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
					      += before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
								   [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
						   * (double)fabs(vesselGraph->vesselSegments[i]->flow) * dt;
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
			     [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
					      -= before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
								   [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
						   * (double)fabs(vesselGraph->vesselSegments[i]->flow) * dt;
		}else{
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
			     [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
					      -= before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
								   [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
						   * (double)fabs(vesselGraph->vesselSegments[i]->flow) * dt;
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
			     [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
					      += before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
								   [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
						   * (double)fabs(vesselGraph->vesselSegments[i]->flow) * dt;
		}*/

		double flow = fabs(vesselGraph->vesselSegments[i]->vesselNodes[0]->pressure - vesselGraph->vesselSegments[i]->vesselNodes[1]->pressure)
				*pow( vesselGraph->vesselSegments[i]->radius, 4)
				/*(double)vesselGraph->vesselSegments[i]->getViscosity()*//DISTANCE(vesselGraph->vesselSegments[i]->vesselNodes[0], vesselGraph->vesselSegments[i]->vesselNodes[1]);

		if( vesselGraph->vesselSegments[i]->vesselNodes[0]->pressure < vesselGraph->vesselSegments[i]->vesselNodes[1]->pressure){
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
			     [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
			      += (before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
							[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
					// -before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
					//		[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
							    						    ) *time_step*flow;
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
			     [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
			      -= (before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
						    [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
					// -before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
					//		[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
									    						    )*time_step*flow;
		}
		if( vesselGraph->vesselSegments[i]->vesselNodes[0]->pressure > vesselGraph->vesselSegments[i]->vesselNodes[1]->pressure){
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
			     [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
			      += (before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
						    [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
					// -before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
					//		[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
									    						   )*time_step*flow;
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
			     [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
			      -= (before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
						   [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
						//	 -before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
						//			[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
											    						   )*time_step*flow;
		}


		/*after[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
		     [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
		      += (//before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
				//	    [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
			     before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
					    [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]])
			   * (double)vesselGraph->vesselSegments[i]->flow * dt;
*/
		if( isnan(after[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
			   		     [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]) ||
			   		  isinf(after[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
			   		  		   		     [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]])
			   		  		   		     ){
			fprintf(stderr, "SHIT! -> %e\n", after[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
			                  		   		     [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]);
			fprintf(stderr, "SHIT! %e\n", before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
			             					    [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]);
			fprintf(stderr, "SHIT! %e\n", before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
			             					    [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]);
			fprintf(stderr, "SHIT! %e\n", vesselGraph->vesselSegments[i]->flow);
			fprintf(stderr, "SHIT! %e\n", flow);
			exit( 0);
		}


		/*after[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
		     [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
		      += (before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
					    [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]])
			     //-before[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
				//	    [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]])
			   * (double)vesselGraph->vesselSegments[i]->flow * dt;
*/
		/*after[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]]
		     [(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]]
		      -= before[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]]
		               [(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]]
		       * vesselGraph->vesselSegments[i]->flow * dt;*/

		if( vesselGraph->vesselSegments[i]->vesselNodes[0]->getType()==ROOT)
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]][(int)vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]] = 1;
		if( vesselGraph->vesselSegments[i]->vesselNodes[1]->getType()==ROOT)
			after[(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]][(int)vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]] = 1;

	}
	temp = before; before = after; after = temp;

	if(it%100==0){
	fp = fopen("oxygen.dat", "w+");
	for (int x = 0; x < DOMAIN_SIZE_X; x++) {
		for (int y = 0; y < DOMAIN_SIZE_Y; y++) {
			fprintf(fp, "%i %i %e %e %e\n", x, y, oxygen[x][y], vegf[x][y], before[x][y]);
		}
	}
	fclose(fp);
	}

	}

	//removeInterTipConnections(vesselGraph, countInterTipConnections);


	exit(0);

	vesselGraph->updatePressure(0, 0, 0);
	vesselGraph->updateShear();

	// SPROUTING
	/*for( int i=0; i<vesselGraph->countVesselNodes; i++){
	 if( vegf[vesselGraph->vesselNodes[i]->position[0]][vesselGraph->vesselNodes[i]->position[1]] > 0.
	 && 1./50. >= RAND01){
	 next = vesselGraph->vesselNodes[i];
	 float maxVEGF = 0.;
	 int candidateX;
	 do{
	 last = next;
	 candidateX = -1;
	 int candidateY = -1;
	 int x = last->position[0];
	 int y = last->position[1];
	 maxVEGF = vegf[x][y];

	 for( int xx=x-1; xx<=x+1; xx+=2)
	 for( int yy=y-1; yy<=y+1; yy+=2)
	 if( xx>=0 && xx<DOMAIN_SIZE && yy>=0 && yy<DOMAIN_SIZE
	 && vegf[xx][yy]>maxVEGF){
	 maxVEGF = vegf[xx][yy];
	 candidateX = xx;
	 candidateY = yy;
	 }
	 fprintf( stderr, "Sprout! maxVEGF=%e\n", maxVEGF);

	 if( candidateX!=-1){
	 next = new VesselNode(candidateX,candidateY,0);
	 vesselGraph->addVesselNode( next);
	 vesselGraph->addVesselSegment( new VesselSegment( last, next));
	 }
	 }while( candidateX!=-1);
	 }
	 }*/


	for (int it = 0; it < 1000; it++) {
		time += time_step;
		updateVesselGraph(oxygen, vegf, vesselGraph, time_step, time);
		updateOxygen(oxygen, oxygen2, vesselGraph, 10000, 1e-5);
		updateVEGF(oxygen, vegf);

		fp = fopen("oxygen.dat", "w+");
		for (int x = 0; x < DOMAIN_SIZE_X; x++) {
			for (int y = 0; y < DOMAIN_SIZE_Y; y++) {
				fprintf(fp, "%i %i %e %e\n", x, y, oxygen[x][y], vegf[x][y]);
			}
		}
		fclose(fp);

		vesselGraph->updatePressure(0, 0, 0);
		vesselGraph->updateFlow();
		vesselGraph->updateShear();

		fp = fopen("vessel.dat", "w+");
		for (int i = 0; i < vesselGraph->countVesselSegments; i++)
			for (int ii = 0; ii < 2; ii++)
				fprintf(
						fp,
						"%i %i %e %e %e %i %lf\n",
						(int) vesselGraph->vesselSegments[i]->vesselNodes[ii]->position[0],
						(int) vesselGraph->vesselSegments[i]->vesselNodes[ii]->position[1],
						vesselGraph->vesselSegments[i]->flow,
						vesselGraph->vesselSegments[i]->shear,
						vesselGraph->vesselSegments[i]->vesselNodes[ii]->pressure,
						vesselGraph->vesselSegments[i]->vesselNodes[ii]->getType(),
						vesselGraph->vesselSegments[i]->radius);
		/*for( int i=0; i<vesselGraph->countVesselNodes; i++)
		 fprintf( fp, "%i %i %e %i\n",
		 vesselGraph->vesselNodes[i]->position[0],
		 vesselGraph->vesselNodes[i]->position[1],
		 vesselGraph->vesselNodes[i]->pressure,
		 vesselGraph->vesselNodes[i]->type);*/

		fclose(fp);
	}

	//vesselGraph->printToPovray( (char*)"test.pov");
	//vesselGraph->updatePressure( 0,0,0);
	//vesselGraph->updatePressure( 0,0,0);

	/*for( int i=0; i<vesselGraph->countVesselNodes; i++)
	 fprintf( stdout, "%i %i %i %lf\n",
	 vesselGraph->vesselNodes[i]->position[0], vesselGraph->vesselNodes[i]->position[1], vesselGraph->vesselNodes[i]->position[2],
	 vesselGraph->vesselNodes[i]->pressure);
	 */

	/*for( int i=0; i<vesselGraph->countVesselSegments; i++)
	 fprintf( stdout, "%lf %lf %e %e %e %e %e\n",
	 0.5*(vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]+vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]),
	 0.5*(vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]+vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]),
	 vesselGraph->vesselSegments[i]->flow*(float)(vesselGraph->vesselSegments[i]->vesselNodes[1]->position[0]-vesselGraph->vesselSegments[i]->vesselNodes[0]->position[0]),
	 vesselGraph->vesselSegments[i]->flow*(float)(vesselGraph->vesselSegments[i]->vesselNodes[1]->position[1]-vesselGraph->vesselSegments[i]->vesselNodes[0]->position[1]),
	 vesselGraph->vesselSegments[i]->flow,
	 vesselGraph->vesselSegments[i]->shear,
	 vesselGraph->vesselSegments[i]->vesselNodes[0]->pressure);
	 */
}
